Skeletal and Dental Development Preserve Evidence of Energetic Stress in the Moose of Isle Royale

Food shortages can leave diagnostic, and in the case of the dentition, irreversible changes in mineralized tissue that persist into historical and fossil records. Consequently, developmental defects of tooth enamel might be used to track ungulate population irruptions or declines in resource availability, but dental tissue’s capacity for preserving historical population density changes has yet to be investigated in wild populations. We test the ability of macroscopic enamel defects, mandible, and metapodial lengths to track changes in the well-known insular moose population of Isle Royale National Park. Our study demonstrates that (1) a moose density threshold exists on the island above for which there is a significant decrease in mandible and metatarsus length and a concomitant increase in enamel hypoplasias; (2) food limitation has a more pronounced effect on male than female skeletal and dental growth; and (3) combined data from tooth enamel hypoplasias and bone lengths reflect the relative density of this ungulate population and should be broadly applicable to other ungulate osteological samples. Developmental defects in dental enamel were among the highest recorded in a wild population, and even during low-density intervals the population density of Isle Royale moose has been high enough to negatively impact skeletal and dental growth, indicating the comparatively poor health of this isolated century-old ecosystem.


INTRODUCTION
A key goal for many paleontologists is characterization of past ecosystems, each aspect of which presents interpretational and taphonomic challenges. The density at which a species existed in its landscape is one such aspect. Recent attempts to characterize evidence of resource limitation in Pleistocene ungulates found evidence of stress following vegetation shifts and probable competition between species for dwindling resources (Barrón-Ortiz et al., 2019). The teeth of Pleistocene carnivorans are also thought to have preserved signals of resource limitation (Van Valkenburgh and Hertel, 1993;Flower and Schreve, 2014) based on studies of tooth wear and fracture in extant carnivorans (Van Valkenburgh, 1988;Van Valkenburgh et al., 2019). However, we know less about how food stress is reflected in preservable elements of ungulate skeletons, and whether it is possible to track the intensity of food stress alongside shifts in ungulate density. We present a neontological study of the potential relationship between events observed in an ungulate population and the development of their mineralized tissue growth with an eye to complementing studies of the recent past. Mammalian teeth and skeletons contain a wealth of information about the life history of individuals and thus can preserve information about populations. Slowed or interrupted growth results from severe metabolic stress, such as might occur in association with an individual's birth, weaning, episodes of disease, and poor nutrition. To test whether skeletal features track changes in population density, we chose ungulates, the prey base of terrestrial large mammalian carnivores for the past 50 million years and of most human societies. Ultimately, we hope to complement existing methods of estimating carnivoran (Van Valkenburgh and Hertel, 1993;Flower and Schreve, 2014) and proboscidean (Fisher, 2001) population dynamics from skeletal material, allowing better characterization of Pleistocene, Holocene, and Anthropocene ecosystems. Such an approach requires characters that can be sampled non-destructively and have a high probability of surviving taphonomic processes, such as teeth. We identified such characters using both tooth and skeletal material from a well-characterized, insular moose (Alces alces) population inhabiting Isle Royale National Park (IRNP), Michigan.
Isle Royale National Park, an approximately 544 km 2 island in Lake Superior, was colonized by moose early in the twentieth century via a rare ice bridge (Murie, 1934). Isolated 90 km off the shore of Ontario, the moose population is unique in that both immigration and emigration are negligible. Since 1948, there has been only a single predator, the gray wolf (Canis lupus), that contributes to moose mortality. The lack of a complete predator guild and steep declines in the wolf and moose populations due to disease have allowed the moose population to fluctuate, punctuated by dramatic population crashes in 1934 and 1996 (Murie, 1934;Vucetich and Peterson, 2020). Over the past 60 years, this isolated predator-prey system has been meticulously documented, facilitating assessment of the effects of climate, population size, disease, and predation on ungulate health and mortality. Throughout this time period, Isle Royale's lowest moose population density has exceeded that of nearby mainland areas, which often fall below 0.03 moose/km 2 [IRNP minimum 0.71; mainland mean 0.55 ± 0.647] (Peterson, 1999;Ripple and Beschta, 2012). When the wolf population has been insufficient to prevent population irruptions, the moose population becomes resource-limited and then rapidly declines rather than stabilizing (Post et al., 2002). Resource limitation resulting from elevated densities of moose has been associated with decreased carcass weight and excessive tooth wear in multiple North American populations (Ferguson et al., 2000), and is argued to be the mechanism driving the rapid dwarfing of the now century-old Isle Royale moose population (Peterson et al., 2011;Silvia et al., 2014).
Several studies have suggested that a variety of skeletal anomalies accompany poor juvenile nutrition in cervids. For example, delayed tooth emergence (Loe et al., 2004;Peterson et al., 2012), shortened metatarsal bones (Geist, 1998;Peterson et al., 2011;Silvia et al., 2014), osteoarthritis (Peterson et al., 2010b), and osteoporosis (Hindelang and Peterson, 1996) have all been linked to malnutrition. Reduced cranial size has been specifically associated with food shortages driven by elevated population density (Hoy et al., 2018). Mandible length is of particular interest because the mandible has a high growth priority over other bones shortly after birth that slows only in the second year of life (Huot, 1988). Prior to the 1996 Isle Royale population crash, shortened mandibles were thought to be linked to elevated population densities (Peterson, 1996) as had been suggested for caribou (Reimers, 1972). However, a study of female moose from another population (Ferguson et al., 2000) found that density affects moose growth rate, not final size.
Developmental defects of tooth enamel have been used to compare population health in two herds of caribou (Wu et al., 2012), but not yet to characterize a population over several decades. The permanent dentition is a prime candidate for preserving evidence of irruptions because tooth tissue does not remodel. The crowns of moose permanent cheek teeth form over about a year, from roughly 3 months before birth to 6 to 8 months post-weaning (Peterson, 1955). Thus, malnutrition can affect the development of tooth enamel during this interval both directly through lack of forage, and indirectly, through an inadequately nourished mother. If foraging juveniles are insufficiently nourished, tooth development will be slowed and/or interrupted, leaving diagnostic features that persist into adulthood. In laboratory and domesticated animals, caloric restriction produces lines of arrested growth in bones, delayed tooth eruption, malocclusion, changes in mandibular shape, and dental enamel hypoplasias (DEHs) (Ablett and McCance, 1969). DEHs are among the most conspicuous, manifesting as a visible localized thinning or absence of enamel (Goodman and Rose, 1990;Guatelli-Steinberg, 2001). DEHs are caused by impaired enamel secretion, a symptom of severe nutritional (Ablett and McCance, 1969;Dobney and Ervynck, 2000), pathological (Suckling et al., 1983;Hillson, 1986;Rose, 1990, 1991) or environmental stress (Mead, 1999;Franz-Odendaal, 2004;Franz-Odendaal et al., 2004;Niven et al., 2004;Byerly, 2007;Wu et al., 2012). The morphology of DEH can vary with the degree of metabolic stress (Witzel et al., 2006), but DEH are easily identified in ungulates with exposed enamel surfaces, such as moose. In wild ungulates, hypoplasias have been attributed to weaning and seasonal stressors in extant and extinct giraffids (Giraffa, Sivatherium) (Franz-Odendaal, 2004;Franz-Odendaal et al., 2004); bison (Bison bison) (Wilson, 1988;Niven et al., 2004;Byerly, 2007;Barrón-Ortiz et al., 2019); horses (Equus spp.) (Barrón-Ortiz et al., 2019); white deer (Odocoileus virginianus) (Davis, 2013); and caribou (Rangifer tarandus groenlandicus) (Wu et al., 2012).
To develop a method broadly applicable to archeological and paleontological samples, we quantified the frequency and severity of growth anomalies (DEH) in the enamel of the cheek teeth. The dental results are compared with the signal preserved in the size of two skeletal elements, mandibles and metapodials, both of which have good preservation potential (O'Connor, 2000). Because this is the first examination of DEH etiology in a wild population, we investigated the interaction between population size relative to resources, predation, and climate. When compared with animals weaned in low-density years, moose that develop during periods of high-population density were predicted to exhibit higher frequencies of anomalies in tissues that form after weaning (enamel, mandibular bone) in addition to shortened stature. If combined dental and skeletal measures reliably track known changes in the moose population, and the development of these characters is conserved in ungulates, then we should be able to use them to identify relative density changes in other populations. We stress that this approach cannot be applied without acknowledging shifts in climate, forage species, disease, predators, and competitors. Even with approximate knowledge of these factors, assessments of nutritional stress using enamel surfaces and bone lengths are primarily a means of gauging population density relative to food availability, and not absolute levels of population density.

MATERIALS AND METHODS
To assess the impact of population density on skeletal and dental characters, we analyzed data on adult moose from IRNP and the adjacent mainland of Ontario (ONT). We studied mandible length, (IRNP n = 940; ONT n = 36), metatarsus length, (IRNP n = 90), and enamel hypoplasias (IRNP n = 137, ONT n = 71). All specimens are listed in Supplementary Table 1. Sample sizes for the mandible, metatarsal, and dental analyses varied in size due to differential availability of specimens. Climate, moose abundance, and predator abundance data are unavailable prior to 1959. Only Isle Royale specimens with complete data were used to generate mathematical models (mandible length n = 940, metatarsus length n = 90, and enamel hypoplasias n = 134).
Skeletal specimens of moose collected from IRNP between 1959 and 2019 are housed in the Michigan Technical Institute Mooseum in Alberta, Michigan. Our study sample spans 1959-2009. The moose of IRNP die of natural causes, primarily predation by wolves but also starvation during severe winters (Vucetich and Peterson, 2004). Samples collected prior to 1980 are primarily wolf kills, whereas deaths after that included a greater number of deaths due to winter severity. After 1980, wolf numbers declined sharply due to disease and did not recover due to inbreeding depression (Vucetich and Peterson, 2004;Robinson et al., 2019). Moose with associated mandibles and metatarsal bones were preferentially sampled. Most IRNP moose have been aged by Isle Royale Wolf Project scientists based on cementum annuli, allowing individuals to be assigned birth years. For 20 IRNP mandible specimens without cementum ages, we estimated a birth year interval based on the specimen's autopsy number (indicative of date of death) and age class at death (e.g., old adult, prime adult). Moose mandible and metatarsus lengths, Isle Royale climate data, and moose density data were downloaded from the IRNP website (Vucetich and Peterson, 2020). For comparison, we also collected dental defect and mandible size data from the closest mainland Ontario populations. The mainland sample (hereafter abbreviated as ONT) spanned the period 1950-1970, and thus experienced similar regional climate conditions as our earliest IRNP low density study period (1959)(1960)(1961)(1962)(1963). ONT specimens are housed at the Royal Ontario Museum in Toronto, Canada. Population density of the ONT population was estimated as between 0.55 and 0.71 individuals per km 2 using Canadian Department of Lands and Forests survey data and maps provided by the Canadian Ministry of Wildlife. The ONT moose population as a whole was subject to gray wolf, black bear (Ursus americanus) predation as well as human hunting, while the IRNP moose population had only gray wolves as predators. To the best of our knowledge, most ONT specimens were killed by humans.

Analysis of Population Density and Morphology
The IRNP moose population has been unstable since the beginning of detailed study in 1959 (Figure 1). Previous work suggests that a density greater than two moose per km 2 can result in stunted bone growth (Ferguson et al., 2000), particularly in the length of the mandible (Peterson, 1996). In addition, the carrying capacity of the island, as determined by browse availability surveys made between 1945 and 1970, was estimated to be 1,000 moose (or 1.83 moose/km 2 ). We used the average of these two density estimates (2, 1.83) as our hypothesized density threshold (1.92) beyond which moose are expected to show evidence of stress. There is evidence that overpopulation and resource limitation can impact moose growth during the year they are born, or during gestation (Palsson, 1952;Solberg et al., 2004) or perhaps even earlier because ungulate population growth has been linked to conditions two years prior to birth (Post and Stenseth, 1998;Patterson and Power, 2002). To use population density as a predictor variable, we needed to establish the year of interest. To do so, we used the maximal information coefficient (MIC) to understand which year metric (abundance at gestation, birth, or two years before birth) was most strongly associated with morphological measures. This let us develop a generalized linear model (ANOVA) with other variables of interest. The MIC assesses both linear and non-linear relationships between variables (Reshef et al., 2011). Analysis was conducted in R 3.3.2 (R Core team) using the package Minerva (Albanese et al., 2013). The analysis was bootstrapped to obtain confidence intervals and p-values for each correlation.
Based on the year with the strongest correlation (birth year, see below), we set the density threshold at 1.92 moose/km and divided IRNP samples into five time periods: two of low density, two of high density, and one transition period (Figure 1 and Table 1). The most recent time period (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) is anomalous due to fluctuating population densities following a 75% reduction in population size. Moose from this time were born under densities ranging from 0.82 to 2.02 moose per km 2 with an average density of 1.86, just below our density threshold. We also believe that the island's ecosystem changed significantly between our early low-density periods and 1997-2008 (see section "Discussion"). We could not subdivide this time period and retain enough samples for meaningful statistical analyses, so we draw most of our conclusions from the first four time periods.

Life History and Tooth Development in Moose
We reconstructed the development of tooth crowns in relation to life history events using published and new radiographs. FIGURE 1 | Moose population density over time. Sampled periods of high density are indicated by shading, as determined using the density threshold 1.92 moose per km 2 (gray line, see text for details). The two horizontal bars (bottom left) indicate the estimated minimum and maximum density of the Ontario (ONT) mainland sample of moose. Note that low density IRNP moose are always at higher density than o bserved for our mainland sample. Life history data were taken from Isle Royale Annual Reports (Vucetich and Peterson, 2020) as well as (Peterson, 1955). We radiographed the mandibles of juveniles of various ages (Supplementary Figure 1) to supplement radiographs from Peterson (1955). The interval of tooth crown development was estimated as beginning in the first month in which radioopaque tissue was visible in the tooth crypt and ending when the full crown was visible to the level of root bifurcation. Tooth development is complete within two years after birth. These ranges were plotted alongside major life history events (Figure 2) to assess which tissues developed near the times of birth and weaning. Differences in male and female moose development that contribute to size differences are also indicated as life history events.

Dental Characters
All surveyed mandibles were radiographed in lateral view to reveal lines of arrested growth and root abnormalities. Two observers (CER, CB) independently surveyed 138 IRNP mandibles as well as 71 ONT mandibles (CB only) for six features that could result from abnormal growth rates in enamel or bone: malocclusion, dental crowding, abnormal wear, tooth fracture, hypomineralization, and DEHs. No ungulate population has yet demonstrated all six simultaneously. Of the six features, only DEHs were observed. This condition is not implicated in chewing impairment or poor body condition as is the case for tooth breakage, heavy wear, and malocclusion when they occur in ungulates (Ericsson and Wallin, 2001;Loe et al., 2006). The widespread pitting, abnormal wear, and characteristic chalky appearance of enamel affected by hyperflourosis (Hillson, 1986;Shupe et al., 1987), another condition that causes enamel pitting in ungulates (Kierdorf et al., 1993(Kierdorf et al., , 1996(Kierdorf et al., , 2000 was absent. Dental enamel hypoplasia were categorized as pits, grooves, missing enamel, and linear enamel hypoplasias following the Fédération Dentaire Internationale Developmental Defects of Enamel Index (Federation Dentarie Internationale, 1982) DDE Index (Supplementary Figure 2). Hypoplasia morphology varies with stress severity (Witzel et al., 2006) and we sought to assess whether hypoplasias manifested differently over time on Isle Royale. Consequently, we assessed hypoplasia data in two ways, first, including all defects, and second, including only the more severe forms. Hypoplasia frequency distribution across groups (five time periods, sex, age) was assessed using bootstrapped chi-squared tests. To analyze differences among time periods, we also used relative risk assessment. Commonly used in epidemiological studies (Zhang and Yu, 1998;Simon, 2001;McNutt, 2003), relative risk assessment estimates occurrence rates of relatively infrequent events and makes pairwise frequency comparisons between groups, such as the moose in each of our six samples (each of the five IRNP time periods plus the mainland).

Mandible and Metatarsus Length
Mandible length was measured as the distance from the most ventral point of the mandibular notch to the most caudal point of the mandibular foramen (Supplementary Figure 2), permitting measurements on specimens with carnivore damage to the angle of the mandible. ONT mandible length data were captured by CB using digital photographs and ImageJ 1 . IRNP adult mandibles more than 4 standard deviation below the mean, 272 mm in length (n = 1), were considered outliers and excluded to facilitate meaningful statistical analysis (median adult length is 347.22 mm; median older adult length is 357.45 mm). 1 http://imagej.nih.gov/ij/

Statistical Analyses
Several factors can confound an assessment of the impact of population density on morphology, such as climate, predator abundance, moose sex, and age (Reimers, 1972). For example, winter severity has a significant effect on moose osteological development (Vucetich and Peterson, 2004;Silvia et al., 2014) and the Isle Royale population's growth rate (Vucetich and Peterson, 2004); thus, it may be expected to play a role in the new characters assessed here. For climate, we used a single index, the North Atlantic Oscillation (NAO) index, an established index of winter severity for Isle Royale (Vucetich and Peterson, 2020) and elsewhere (Stenseth et al., 2003). On Isle Royale, winters tend to be more severe when the NAO trends negative. Data on snow depth, the NAO and species abundance were derived from the Isle Royale Wolf Project website (Vucetich and Peterson, 2020). Moose sex and age can be expected to interact, because adult male moose continue to grow for several years past sexual maturity while females do not (Garel et al., 2006). We analyzed DEH counts and bone lengths as pooled samples as well as subdivided by sex. We first obtained a variance inflation factor to identity variables that were collinear using constrained correspondence analysis within the R package rioja (Juggins, 2015;R Core Team, 2019). To explore the relative impact of the remaining factors on our morphological characters, we assessed general linear models (GLMs) that included log-transformed data of moose age, a climate proxy, moose density, and wolf density and analyzed them using the RRPP package in R (Collyer and Adams, 2018). Model performance was evaluated based on log likelihood scores and the relationship and significance of terms in a multivariate model were evaluated using their coefficients and multiple R squared values. Because our bone length data do not follow assumptions of normality (mandible length, Shapiro-Wilkes W = 0.79, p < 2.2e−16, metatarsus length W = 0.89, p-Value = 9.86e−07), we used a bootstrap in all analyses and chose tests robust to such data. To assess the effect of population density on moose mandible length between sexes and age groups, we used a three-way ANOVA. To further analyze significant interactions between pairs of terms, such as sex and population density, we used two-way ANOVAs.
All code used in these analyses, as well as dental defect data are available by request to CB. Measurement data for ONT moose can be obtained from datadryad.org. IRNP mandible and metatarsus length data are property of the Isle Royale Wolf Project.

Generalized Linear Models
Variance inflation factor analysis showed that the two forms of climate data were strongly colinear, with snow depth contributing more to this issue than the NAO. Snow depth data are missing for years prior to 1974, thus snow depth was eliminated, leaving population of predators, population of prey, age in years, and the NAO in the models. These models typically explained a small amount of the variance in bone length data ( Table 2 and  Supplementary Table 2). No models of DEH estimated as the number of DEH per individual were significant and the majority had goodness-of-fit estimates less than 0.02. When evaluating all moose together, the best models included one or more of three variables, age, wolf abundance, or moose abundance, although all terms in the best model were not always significant (Supplementary Table 2). When the results were divided by sex, the best-performing models for male jaw and metapodial lengths included age and climate in addition to predator and prey density; models for female bone lengths rarely included climate as significant. Age was significant in most mandible and metatarsus models indicating that there is some growth beyond sexual maturity in both males and females. It is notable that age also appeared in some of the best models predicting DEH even though moose teeth develop at the same age in all moose. Overall, the influence of predator abundance, moose abundance and age were clear, although variance in the data was very high and the number of DEH could not be modeled. Subsequent analyses (ANOVA, chi-squared tests) examined moose within each period, with sexes and younger adults separated. DEH were analyzed by the proportion of the population with the character rather than by number of DEH.
Some variables available to us from Isle Royale would be impossible to measure precisely from the fossil record such as yearly climate change, annual counts of DEH per individual in a population, and precise predator numbers. To obtain predator abundance relative to resources, methods based on tooth fracture and wear might be used to estimate relative abundance, but cannot provide absolute estimates like the wolf numbers in the models included here (Van Valkenburgh and Hertel, 1993;Flower and Schreve, 2014). It is likely that low predator numbers are the ultimate cause of resource limitation on the island; unsustainable moose population size is the proximate cause of bottom-up limitation. Because our intent was to recover possible proxies for energetic stress in ungulates that can be read from collections of time-averaged fossil specimens, we focused our remaining analyses on the proportion of moose in each time period with one or more DEH and their mean bone lengths.  -5.943, 5.667], p < 2e−16). They also exhibited elevated frequencies of dental anomalies (DEH), relative to moose born during the two low density periods (Figures 4, 5 and Table 1). Mandible length of IRNP moose varied only slightly over the first three time periods. It then declined sharply in period 4, the second of the high-density periods. Metatarsal length also shows its steepest decline between periods 3 and 4 (Figure 4). The differences in size and DEH frequency over time suggest that overpopulation on IRNP had a negative impact on both growth of teeth and bones. This is further supported by the comparison of mainland and IRNP moose. The frequency of DEH in IRNP moose ranges from 33% (period 1) to 81% (period 2) with a mean of 62.45% as opposed to 36% in the mainland moose sample. After period 1, the frequency of DEH in IRNP moose was 40-200% greater than that observed in the ONT sample. Moose body size, measured by mandible length, was reduced on Isle Royale in the earliest period of study when compared to contemporaneous mainland moose. Although metatarsal length data was not available for the mainland sample, the moose residing on IRNP in period 1 had significantly shorter mandibles than the Ontario sample (effect size = 20.53 mm, 95% CI [7.1155, -10.51], p = 1e−04). However, the sample size of measurable Ontario mandibles was small (n = 36) and skewed toward males, and the sex ratio of the IRNP period 1 sample is skewed toward females (1.68:1), which accentuates the size difference. Individuals with shorter bone lengths tended to have more DEHs, and these were not evenly distributed across teeth and cusps (X 2 = 473.37, df = 5, n = 418, p < 0.001). The premolar teeth are nearly free of hypoplasias whereas certain molar teeth, such as the m1 and m2, are particularly prone to them (Figure 5). Of all moose exhibiting enamel defects, 36.8% (77/209) developed a linear defect on the buccal surface of the m2 approximately 3 mm above the cemento-enamel junction while 22.9% (48/209) developed a linear defect on the buccal surface of m1. The  This implies that stressors other than post-weaning nutrition contribute to the total hypoplasia count. Likely candidates are some earlier aspect of life history, such as birth or weaning (see Figure 2).
Dental enamel hypoplasias were recorded in all IRNP periods, but the frequency distribution is significantly uneven over time (X 2 = 18.521, df = 5, n = 209, p < 0.001). DEH frequency peaks in period 2 declines in period 3, and then rises slowly across periods 4 and 5 (Figure 4). This pattern is driven by male moose; the proportion of females with enamel hypoplasias did not change significantly across periods (males X 2 = 9.5065, df = 4, n = 45; p-Value = 0.04961; females X 2 = 4.9811, df = 4, n = 33, p-Value = 0.2892). Risk of DEH in IRNP moose was also elevated during both high-density periods relative to lowdensity periods, and relative to the ONT sample ( Table 3). If premolars are formed entirely after weaning, then the relative risk of forming a defect increases once a moose stops nursing and forages. IRNP moose from high-density period 2 had the greatest risk of forming DEH, particularly post-weaning hypoplasias that can be directly linked to insufficient forage. However, the risk of more severe pit and missing enamel DEH morphologies was not significantly elevated in any period ( Table 3), demonstrating that while presence of DEH was associated with population density, DEH morphology was not.

Effect of Age and Sex of Mandible and Metatarsal Lengths
As revealed by our models, age, sex, population density, and the density of predators, all had a significant effect on the length of the mandible (Table 4). Pairwise investigations of sex, population density, and bone length showed that high population density had a significant effect on median male mandible length (effect size = -5.28mm, 95% CI [-3.08, -9.33], p < 0.001) but not FIGURE 5 | Proportion of moose exhibiting one or more hypoplasias per tooth and time period or locality. Teeth are displayed in approximate order of development and eruption. Gray and black bars indicate high-density periods of study on Isle Royale. Ontario moose represent a mainland sample contemporaneous to period 1. The relative risk of pit and missing enamel DEH reflects the occurrence of more severe and longer-lasting stress than the relative risk of all DEH types. IRNP, Isle Royale National Park; ONT, closest mainland Ontario moose.  Figure 3); however, the size difference between the sexes is halved under high-density conditions (median difference at low density = 10.29 mm; high density = 5.33 mm). Similarly, moose sex, population density, and the interaction between sex and density were significant contributors to adult metatarsus length. Two-way analysis of metatarsus length, sex, and density showed that population density, but not sex, had a significant association with metatarsal length on its own. The fact that sex alone was not significantly associated with metatarsal length is surprising, especially in view of the overall size difference between males and females and the significant association with mandible lengths (Figure 4 and Table 4) (mean difference 6.157 cm, [-4.84, 16.521] p = 0.5118).
Although final metatarsus length is strongly influenced by in utero nutrition (Palsson, 1952), adult metatarsus lengths were not significantly shorter in individuals with greater total numbers of hypoplasias or greater numbers of hypoplasias incurred prior to weaning.

DISCUSSION
Our findings strengthen the proposed links between population size, food limitation, and growth interruptions in ungulate teeth and skeletons. Using a wild population of moose, we found strong associations between predator decline, population irruptions, resource limitation, and DEH. Sustained low predation rates and elevated moose population size were the most important contributors to reduced stature and DEH formation in the population as a whole and particularly in male moose. Both dental enamel and bone lengths reflected long-term population trends in the Isle Royale moose population. Our analyses confirm that a moose density threshold exists on IRNP between 1.8 and 2 moose per km 2 , above which there is a significant decrease in mandible and metatarsus length and a concomitant increase in enamel hypoplasias. Densitydependent decreases in moose carcass weight, and presumably skeletal size, have been previously demonstrated (Ferguson et al., 2000;Solberg et al., 2004), as well as increases in incisor tooth breakage (Clough et al., 2010). Here we confirm that signals of overpopulation may be preserved in both bones and teeth, raising the possibility of estimating relative population size from skeletal materials. The frequency of post-weaning enamel defects appears to reflect food-limitation. When sexes are pooled for the IRNP moose, the relative risk of post-weaning enamel defects is significantly higher under high-density conditions. There is a conspicuous near-absence of post-weaning hypoplasias in both low density, pre-crash (period 1) Isle Royale and mainland Ontario moose. Bone lengths are always significantly reduced on the island, but the effect is more pronounced when population density is high as a result of low or effectively absent wolf predation.
The number of DEH per individual did not have a linear relationship with moose density or wolf density. This is not unexpected because a short time of stress could affect as few as one or as many as four teeth that are developing during the same window (see Figure 2); consequently, a greater number of DEH in one individual does not necessarily indicate a higher level of stress across the population. Because we are interested in an estimate of population rather than individual health, the proportion of the population exhibiting one or more hypoplasias is a more informative measure, and this number varied strikingly between periods. After 1968, the frequency of DEH in IRNP moose was much greater (40-200%) than that observed for mainland Ontario moose. Because our IRNP moose sample is dominated by wolf-killed moose, and because wolves select individuals that are easier to kill, it is possible that our IRNP sample is dominated by relatively unhealthy individuals. If the presence of DEH is indicative of poor health in adults, then the very high frequency of DEH we observed in IRNP moose relative to Ontario moose might reflect the difference between wolf-selected as opposed to human-selected moose. However, the presence of DEH is unlikely to be indicative of adult health given that they form in the first two years of life, and these small features are not associated with chewing impairment or poor adult body condition in ungulates. It is much more likely that the high frequency of DEH in IRNP moose record nutritional or other metabolic stresses early in life without contributing to adult mortality. Moreover, we see clear signals of metabolic stress within some time intervals that reverse rapidly (i.e., period to period) when moose density declines and competition for resources lessens. Period-by-period examination of DEH rates shows them to be more prevalent when wolf numbers are low, including during period 2 when the wolves' ability to take prey was not obviously compromised by disease or inbreeding depression. In our GLM analyses, predation was the most consistent and important predictor of metatarsus length; unlike enamel furrows, leg length affects locomotor ability. Individuals of shorter stature will experience lasting impairments, possibly explaining the greater importance of wolf predation in models that predict metatarsus length as opposed to jaw length.
Comparison of moose on and off the island provided interesting insights into the health of the Isle Royale ecosystem over the past 70 years. Our combined study of tooth and bone suggests that even when population densities were at their lowest on the island, IRNP moose experienced significantly more nutritional stress than their contemporaneous mainland counterparts, where the density was five to ten times lower. The earliest-studied Isle Royale moose had 1.18 times the chance of incurring a hypoplasia and 2.98 times the chance of incurring a post-weaning hypoplasia relative to their contemporaries on the mainland. Period 1 IRNP moose also had significantly shorter mandibles than the mainland moose although this might reflect the larger proportion of females in our IRNP period 1 sample. However, our results corroborate previous observations that moose began rapidly dwarfing within decades of their early 1900's arrival on IRNP (Silvia et al., 2014) and are now smaller in size than mainland Michigan moose (Peterson et al., 2011). The rapid dwarfing and elevated DEH suggest that the IRNP moose were suffering from limited forage on the island relative to the mainland soon after their arrival.
Once on the island, moose exerted heavy browsing pressure that altered forest composition through the suppression of woody plants (Murie, 1934;Risenhoover and Maass, 1987;McInnes et al., 1992;Rotter and Rebertus, 2015) and extirpated important browse species (Murie, 1934;Brandner et al., 1990;Vucetich and Peterson, 2004). It takes more time for these species to recover than it takes for the moose population to rebound (McLaren and Peterson, 1999). In particular, local shortages of the winter browse species balsam fir (Abies balsamea) (Brandner et al., 1990;McLaren and Peterson, 1999) may be linked to the prevalence of hypoplasias on late forming teeth in periods 2-5 (Figure 5), regardless of moose population density. Our results underscore the need for sustained periods of low moose population density to support both the moose and the vegetation on which they depend. This will only be feasible if the wolf population grows, but it may become a reality following the National Parks Service's recent decision to repopulate the island with wolves (Mlot, 2018) and given the lower than expected moose population growth under atypically warm climate conditions (Vucetich and Peterson, 2020). The continued decline of the island ecosystem's browse availability and the health of the wolf population (Peterson et al., 2010a;Robinson et al., 2019) likely contributed to the poor fit of our mandible and metatarsus models.
Signals preserved in bone and enamel demonstrate a clear effect of elevated IRNP moose density and low predation at the population level and between the sexes. While the overall frequency of DEH was elevated under high-density conditions, the risk of more severe pit and missing enamel DEH morphologies was not, indicating that overall DEH frequency is more informative than relative numbers of different types of DEH. DEH differed significantly in frequency between both sexes and population density regimes, with density-related changes driven almost entirely by males. Whereas Davis (2013) found no sex differences in DEH frequency in white-tailed deer, we found that male IRNP moose had significantly elevated frequencies of DEH during high-density periods. In contrast, female DEH frequency in IRNP moose varied much less over time. Moose generation time is short (Geist, 1963;Van Ballenberghe and Ballard, 1994;Gaillard, 2007), and male reproductive success is driven by body size (Garel et al., 2006). Consequently, male dental growth seems more affected by conditions early in life than females, due to the combined tendencies of males to grow faster and for a longer period of time (Kruuk et al., 1999;Solberg et al., 2004;Garel et al., 2006). In moose, it seems that the sexes differ in their sensitivity to metabolic stress as early as birth, the time at which the first permanent molar, the earliestforming feature in our study, develops. The presence of age in DEH models is unexpected, but may reflect the disappearance of DEH as tooth crowns wear down. We selected very old individuals with intact crowns and enough enamel to evaluate the most commonly affected areas of the teeth (see below). However, our results indicate that very old age or extreme tooth wear is likely to cause DEH to be underestimated. Moose are more likely to reach very old age (>10 years) when predator abundance is low.
Across all five IRNP time periods and both the island and mainland populations, enamel hypoplasias appeared in constrained locations of the tooth crown that develop around the time of birth (m1) and weaning (m2). During high-density periods, there was a significantly greater risk of developing an enamel hypoplasia on tooth crowns formed both prior to and during weaning. This suggests that the development of DEH during weaning is most likely to occur in the presence of population-level resource limitation, although they can also occur in its absence. This is similar to the phenomenon reported by Wu et al. (2012) in which DEH associated putatively with birth were more common in a caribou population experiencing a population decline than non-declining populations. Despite the effect of a localized thinning of tooth enamel, DEH on molar teeth were not associated with broken teeth or abnormal wear, and thus we have no reason to suspect that DEH caused individuals to forage less efficiently; as noted above, they are a signal of metabolic stress in the past and not a contributor to poor condition in adulthood.
The peak hypoplasia frequencies observed on Isle Royale (81%) are among the highest reported in wild populations of ungulates, with severe impairment of enamel formation visible on many individuals. A nutritionally stressed population of white-tailed deer exhibited defects in 27% of individuals (Davis, 2013), while a population of caribou in decline had a frequency of 25.1% (Wu et al., 2012). Although a declining population of Alaskan moose exhibited higher levels of enamel defects (93%) (Stimmelmayr and Persons, 2006), the defects were not localized lesions such as the DEH described in IRNP moose and the above populations, but rather structural failures affecting most of the tooth crown that were putatively linked to local geochemistry (Clough et al., 2010;MacKenzie et al., 2011). Given that IRNP moose do exhibit reduced heterozygosity (Wilson et al., 2003), it is possible, but unlikely, that their extremely high frequencies of DEH result from inbreeding. To date, DEH frequency has not been shown to increase in inbred populations; furthermore, known genetic enamel disorders are characterized by thinning and pitting of all enamel surfaces, abnormal tooth shape and malocclusion in humans (Hillson, 1986;Wright et al., 1993) and mouse models (Gibson et al., 2001;Caterina et al., 2002;Poulter et al., 2014). In our IRNP moose samples, DEH appeared on otherwise normal teeth with no signs of abnormal wear due to malocclusion or other factors.
Unlike moose teeth, which complete their growth within a year or two at most, bones have the potential to capture the effects of food limitation over a longer period of time. Moose mandible and metatarsus lengths varied with population density as expected from existing life-history data; higher population density reduces food availability and consequently moose stature. Like studies of moose populations elsewhere (Ferguson et al., 2000;Garel et al., 2006), we found that high population density has a much greater effect on male than female size. Studying exclusively females, Ferguson et al. (2000) documented shortened mandibles in growing Canadian moose, but final adult mandible sizes did not differ under high-density conditions. We were able to quantify mandible lengths in both sexes and found a significant reduction in males but not females consistent with Ferguson et al. (2000). The greater impact of population density on males may again be attributed to their longer period of growth. Just as we observed from the tooth enamel, male skeletal growth slows while they are stressed, whereas females may compensate through other means such as delayed breeding (Reimers, 1972;Sand, 1996;Keech et al., 2000;Testa, 2004;Garel et al., 2006).
Our study of IRNP moose demonstrated population-level changes in morphology but also found discrepancies between tooth enamel and bone length data. Mandible and metatarsus length tracked population density more closely than did DEH frequency (Figures 4, 5). Within high-density periods, the maximum population density is matched by minimum average bone lengths. Unexpectedly, moose from period 2 had the most frequent DEH and the greatest risk of incurring a DEH, suggesting an extended period of stress, but not the greatest reduction in bone lengths. Moose born immediately post-crash in period 5 had similar DEH to moose born under peak density, but recovered mandible lengths. These contradictions may result from the longer growth period of bones in contrast to teeth. Because mandibles and metatarsi continue to grow for several years after birth, a reduction in growth in earlier years might be reversed in later years if food becomes plentiful. If this approach is to be applied to other systems, it is important to note that compensatory growth (or catch up) later in life has been demonstrated in cervids and other ungulate species (Luke et al., 1981;Suttie et al., 1983;Watkins, 1991;Freetly et al., 2000;Betterly, 2011) and could obscure the population density signal found with our approach.
Unlike bone length, which varies significantly between individuals and even cohorts of individuals on a similar nutritional plane (J. Vucetich, pers. comm.), teeth complete their growth at specific times and are relatively constrained in size. Consequently, tooth enamel and DEH frequency should provide a more reliable record of metabolic stresses early in life than either mandible or metatarsal lengths across ungulate species as long as wear has not removed most of the tooth crown. These lesions possibly correspond to stressors other than food limitation-a phenomenon we observed in period 2, a time characterized by severe winters and peak DEH but not population size. Stressors during this period could have included energetic expenses not captured by the NAO index. Examples are winter ectoparasites (DeIgiudice et al., 1997) or locomotor difficulties resulting from a dense winter snow crust. Parasite outbreaks that coincided with times of high density (particularly in 1987, 1989, and 1996) would have added to the number of individuals with hypoplasias in period 4 as well. There will always be factors, such as disease and parasites, that are unlikely to be apparent in the fossil record, but we can observe relative changes and evaluate multiple hypotheses regarding resource shortages by including both bones and teeth, climate data, and estimates of health in other taxa.

Application to the Fossil Record: Another Tool in the Kit
The present study is a first-step neontological study intended to discover the relationship between relative population density and skeletal morphology of elements often found in the fossil record. Given that we conducted this study with a goal of applying our findings to ungulate populations throughout the historical and paleontological record, the likelihood of preservation of characters is just as important as their ability to record episodic stress. In the fossil record, dental hypoplasia frequency, particularly during post-weaning, is likely to preserve a record of resource limitation or other metabolic stress that is also recorded in mandible or metatarsal length. In addition, teeth do not continue to grow or remodel and are often better preserved than bones. Therefore, skeletal measures should be considered in conjunction with dental development, particularly in the context of prehistorical applications. We stress that the measures presented here provide only relative estimates of population density. Nonetheless, if our study of Isle Royale is viewed on the coarser time scale of a century, the evidence of periods of overpopulation and very low predator to prey ratios are preserved consistently, while shorter intervals of difficult winters or disease would not be apparent. It is probable that vegetation changes would preserve in paleoecosystems that experienced conditions similar to Isle Royale because of the decline or extirpation of some browse species. These localized changes could also be attributed to an excess of moose herbivory from overpopulation rather than climate change. This gives us confidence that this same ecosystem dynamic could be recovered from the past given a chronological sequence of preserved skeletal elements.
Our study of IRNP moose has established a likely correlation between dental defects and overpopulation that adds an exciting new tool to our kit of morphological indicators of population stress that can be applied to fossil ungulates. The ability to draw evidence from multiple fossilized tissues when estimating relative population size may be particularly useful in reconstructing the terminal Pleistocene when predator guilds collapsed. Assuming that the frequency of dental enamel defects reflects relative population health in fossil ungulates, it might be possible to track trends in Pleistocene populations of ungulates over hundreds or thousands of years, as they approach their extinction, as well as comparing DEH frequency among survivors before and after the removal of megafaunal predators and competitors. Analysis of dental development in conjunction with climate proxies can indicate whether populations of a species might have declined due to bottom-up (e.g., climate) limitations, as opposed to top-down (e.g., predation) pressure (see Fisher, 2001;Barrón-Ortiz et al., 2019). After the late Pleistocene extinction, surviving ungulates may have experienced a release from predation pressure and interspecific competition, allowing their numbers to increase beyond carrying capacity. In such a scenario, irruptions might be manifest in increased DEH frequency and shorter mandibles and metatarsi. When combined with data on tooth wear and fracture in co-existing predator species (Ripple and Van Valkenburgh, 2010;Van Valkenburgh et al., 2015), it has the potential to provide new insights into the predator-prey dynamics of ancient ecosystems.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the DataDryad repository (https://datadryad.org/stash/share/ lGt737JG9KrGyDJ4ctMrGD1DSCOtyRz1Fj4GLtJt2MY).

AUTHOR CONTRIBUTIONS
BV, WR, and CB conceived the project. BV, CB, and CR designed the research and collected the data. CB executed the statistical analyses and wrote the first draft of the manuscript. BV, CR, and WR contributed to the writing of the manuscript. All authors gave final approval for publication.

FUNDING
This work was supported by the National Science Foundation EAGER Award Number 1237928 to BV.