75 TWIG PRODUCTION PATTERNS AMONG MOOSE FORAGE SPECIES AND IMPLICATIONS FOR FOREST MANAGEMENT Matt T. Petz Giguere1, William J. Severud2* Kim Teager 3†, Tiffany M. Wolf2‡,3 Seth A. Moore4‡ 1Forestry Department, Gichi Onigaming Anishinaabeg (Grand Portage Band of Lake Superior Chippewa), Grand Portage, Minnesota 55605, USA; 2Veterinary Population Medicine, University of Minnesota, St. Paul, Minnesota 55108, USA; 3Natural Resources Management, Lakehead University, Thunder Bay, Ontario, P7B 5E1, Canada; 4Biology Department, Gichi Onigaming Anishinaabeg (Grand Portage Band of Lake Superior Chippewa), Grand Portage, Minnesota 55605, USA Correspondence author: Matt T. Petz Giguere, Forestry Department, Gichi Onigaming Anishinaabeg (Grand Portage Band of Lake Superior Chippewa), Grand Portage, Minnesota 55605, USA. Email: mtyler@grandportage.com ABSTRACT: A declining moose (Alces alces) population in northeastern Minnesota is a serious con- cern of the Gichi Onigaming Anishinaabeg (Grand Portage Band of Lake Superior Chippewa). A better understanding of the variation in moose forage production at the plant level among tree and shrub species and over time may be important to moose recovery in areas where nutrition may be limiting, while managers need practical methods to monitor forage density. Using data from an exten- sive 2019 moose browse survey at Gichi Onigaming (Grand Portage Indian Reservation, MN) and Minong (Isle Royale, MI), we fit models of twig production (number of twigs per plant) within moose reach as a function of species, canopy cover, and stem height. We validated the best fit model against an independent data set and used the model to simulate stand-level forage density over time for 3 tree species preferred by moose. Twig production varied non-linearly with stem height. An increasing then decreasing unimodal curve with height fit better than allometric models or species means. Peak twig production, height at peak production, and rate of production decline with height varied among spe- cies. Paper birch (Betula papyrifera) and balsam fir (Abies balsamea) generally had the greatest peak twig production, greatest heights at peak production, and lowest rates of decline with height. Among trees, quaking aspen (Populus tremuloides) generally had lower peak twig production, lower height at peak production, and greater rates of decline with height than most other trees and some tall shrubs. Peak twig production was greater under open canopy than under closed canopy for seven species. Model validation indicated that predictions were well correlated with observations and outperformed alternative models but had consistent over-prediction bias. Simulations of regenerating aspen, paper birch, and red maple (Acer rubrum) forests indicated that whereas aspen produced ~ 1.4-2.2 times more peak forage biomass than paper birch or red maple, paper birch and red maple produced usable forage densities for ~2-17 years longer than aspen. Finally, we provide a procedure to use the regres- sion equations to estimate moose forage density from common forest regeneration survey data. Although the equations are suitable for monitoring, they should be used cautiously in high accuracy applications. Our findings suggest that diverse mixes of deciduous trees and shrubs resulting from post-harvest treatments likely provide abundant moose forage for longer durations than do nearly pure aspen stands often resulting from clearcutting alone in northeastern Minnesota. ALCES VOL. 60: 75–108 (2024) Key Words: Alces alces, forage equation, indigenous, Minnesota, moose, twig production *Present address: Department of Natural Resource Management, South Dakota State University, Brookings, SD 57007, USA †Present address: Lake Superior National Marine Conservation Area, Nipigon, Ontario P0T 2J0, Canada ‡Co-senior authors mailto:mtyler@grandportage.com TWIG PRODUCTION PATTERNS ALCES VOL. 60, 2024 76 Moose (Alces alces; Mooz in Anishina- abemowin, the Anishinaabeg language) are a primary subsistence species used by the Gichi Onigaming Anishinaabeg (Grand Portage Band of Lake Superior Chippewa), a native community in northeastern Minnesota who exercise rights to hunt, fish, and gather across their traditional homeland (Stults et al. 2016). Maintaining subsistence species into perpetu- ity through the Anishinaabeg practice of sev- enth-generation planning (Vukelich 2023) and scientific inquiry are of high priority to the Tribal nation (Moore et al. 2024). As one of several strategies to maintain or increase moose numbers, the Gichi Onigaming Forestry Department conducts forest management proj- ects at Gichi Onigaming intended to increase the abundance of woody forage preferred by moose. The data and research presented in this article were gathered with and for the Gichi Onigaming Anishinaabeg. Moose populations in Minnesota have been declining (>50% decline since the 1990’s) (Lenarz 2007, Severud et al. 2022), and population declines have impacted the Gichi Onigaming community’s traditional way of life. To address these concerns, the Gichi Onigaming Anishinaabeg began a long-term moose research program in 2010 and have found that the decline has resulted from low survival rates of adults and calves (Wolf et al. 2021, Severud et al. 2022, Wehr et al. 2024, Garwood et al. 2025). Research in Minnesota during 2003–2022 by the Minnesota Department of Natural Resources (DNR) and others have indicated that para- sites, disease, predation, and winter nutri- tional restriction are contributing to moose population declines, resulting in high mor- tality of adults and poor recruitment of calves (Severud et al. 2015, 2019, 2022; Wünschmann et al. 2015; DelGiudice et al. 2019; Carstensen et al. 2018). Management to improve forage condi- tions may alleviate the negative effects of nutritional limitation on moose population dynamics in this region. Recent studies in other regions suggest that forest manage- ment decisions positively or negatively affect the abundance of preferred forage (Johnson and Rea 2024) and that the result- ing “nutritional landscape” may contribute to moose population growth or decline (Schrempp et al. 2019; Peterson et al. 2020, 2022). Studies have found positive associa- tions between preferred forage availability and calving success (Hayes et al. 2022) and between moose calf weights and a diverse diet of deciduous species (Felton et al. 2020). In Adirondack Park, New York, USA, moose densities were positively associated the abundance of summer forage with high digestible protein (Peterson et al. 2022), which was often a result of forest harvest. Also in Adirondack Park, local browsing intensity was associated positively with the density of preferred forage plant species and negatively with that of avoided species (Peterson et al. 2020). Together, these stud- ies indicate that the abundance of diverse, high-quality preferred forage can influence moose at the population level. In northeastern Minnesota, forest man- agement practices are often directed at com- mercially valuable quaking aspen (Populus tremuloides; Anishinaabe language: Azaadi; De Pellegrin Llorente et al. 2020), which is preferred by moose (Portinga and Moen 2015). Aspen in northeastern Minnesota is commonly managed in even-aged stands harvested by clearcut with reserves (Perala 1977, David et al. 2001). Aspen typically regenerates after clearcutting through high densities of root suckers that grow rapidly (Perala 1977, David et al. 2001) and initially provide abundant moose forage. However, little is known about summer protein content in quaking aspen leaves relative to other for- age species (but see Renecker and Hudson (1988)), and aspen forage quality is further ALCES VOL. 60, 2024 TWIG PRODUCTION PATTERNS 77 complicated by seasonal variation in tannins (Lindroth 2001) that reduce protein avail- ability (Robbins et al. 1987). Similarly, aspen’s rapid height growth allows it to out- compete other forage species (Zenner and Puettmann 2008), reducing diversity while also potentially reducing its own accessibil- ity over time. Recent research also suggests moose may prefer burned areas over har- vested areas (Schrage 2023, Mumma et al. 2024), although the mechanisms are not clear. Overall, it is possible that modifica- tions to forest management might produce more diverse, higher quality, or longer last- ing forage for moose. To improve the nutritional landscape for moose, managers need a better understanding of forage dynamics as well as rapid and prac- tical methods to estimate forage abundance and quality. Ideally, forage quality and quan- tity could be incorporated into existing moni- toring programs, such as post-harvest forest regeneration surveys. Such methods may need additional protocols to accurately esti- mate browse abundance that is available to moose. Moose typically eat twig ends during the dormant season and stripped leaves during summer (Renecker and Schwartz 2007). Twigs and leaves between 0.5 m to 2.5-3.0 m above the ground are used most heavily, pre- sumably because they are easiest to reach (Renecker and Schwartz 2007). Although moose sometimes strip bark or push over large tree saplings to access twigs above 3 m, this is usually associated with exceptional conditions such as low food availability or difficult winters (Peek et al. 1976; Renecker and Schwartz 2007). The only study to estimate browse pro- duction by individual woody plants for northeastern Minnesota (Grigal and Ohmann 1977) estimated leaf and twig biomass using measurements of plant diameter at 15 cm above ground for shrubs and trees. However, the power law equations used for these estimates do not account for whether twigs are within reach of moose, and many have poor coefficients of determination. For example, the Grigal and Ohmann (1977) equations estimate that a mature aspen tree 25 cm (~10 inches) in diameter at the base would have ~2,800 grams of current annual twig biomass. However, it is reasonable to expect that nearly all of that biomass would be well out of the reach of moose and thus unavailable. Increasing and then declining browse availability with the age of forage plants is consistent with multiple research studies showing that moose use of harvested or burned areas increases after the distur- bance and then declines after about 10-20 years (Loranger et al. 1991, Fisher and Wilkinson 2005). Mechanisms behind this pattern may include changes in species com- position as well as changes in plant height that create a disconnect between total bio- mass and accessible biomass. Instead of Grigal and Ohmann’s allometric equation, the relationship between browse availability and plant size might be better modeled by a unimodal function that rises to a peak and then declines again to low values, which would account for reduced browse availabil- ity to moose beyond an upper limit in plant height. For example, similar browse biomass equations with unimodal properties have been recently developed in another ecologi- cal system within Adirondack Park (Peterson et al. 2022). The Grigal and Ohmann (1977) or Peterson et al. (2022) equations would also be difficult for foresters to apply to common field data because field foresters do not often measure plant basal diameter. Instead, foresters often collect tree regenera- tion and shrub competition data in terms of density by species and height (e.g., Minnesota DNR Forestry regeneration mon- itoring procedures [Minnesota DNR Forestry 2016], Wisconsin DNR forest regeneration metric [Wisconsin DNR 2021]). TWIG PRODUCTION PATTERNS ALCES VOL. 60, 2024 78 Multiple studies have quantified species–specific characteristics of twigs eaten or available to be eaten by moose. Portinga and Moen (2015) provided tables of mean oven dry mass and diameter of twigs used by moose in northeastern Minnesota for 12 spe- cies in winter and 11 species in summer under open and closed canopies. Peek et al. (1976) provided estimates of mean browsing diame- ter and mean oven dry mass per twig of 10 species of forage eaten in fall and winter in northeastern Minnesota. Risenhoover (1987) provided estimates of mean browsing diame- ter and mean oven dry mass per twig of 17 species of forage eaten in winter in Isle Royale National Park (Minong). Given the availabil- ity of individual twig biomass data, determin- ing the relationship between woody plant characteristics and the number of available twigs produced per plant may allow for more efficient models of biomass estimation per plant, and might also integrate seasonal and species corrections. Such an approach would not only account for biomass available to moose, but also meet the need for efficient data collection by foresters and managers. Prior research also suggests that plant- level browse production may vary between species and with differences in canopy clo- sure or browsing pressure. The variability among the coefficients of the Grigal and Ohman equations suggests production differ- ences between species. De Jager and others working at Minong and in Sweden found that tree species differed in their response to browsing; some species increased twig pro- duction in response to moderate moose browsing whereas other species decreased production (De Jager 2008, 2009; De Jager and Pastor 2008). Portinga and Moen (2015) found differences in individual twig biomass between species and between open vs. closed canopy conditions within species. It is there- fore possible that canopy closure may simi- larly affect the number of twigs produced. To better understand the relationship between woody plant characteristics and the production of twigs available to moose, we developed a series of models using data from extensive moose browse vegetation surveys at Gichi Onigaming and Minong. We then verified the models against an independent data set and conducted a simple simulation of the forest management implications of our findings. Finally, we offer a brief worked example of how our models can be applied by field staff to rapidly estimate forage availabil- ity in the field. We hypothesized that: H1. Plant-level twig production within the acces- sible browse zone by trees and tall shrubs has a complex non-linear relationship with plant height, first increasing with height as the plant grows more branches, but then decreasing with height as lower branches senesce and new twigs are out of reach. H2. Plant-level twig production curves vary among species. H3. Plant-level twig production curves vary with forest overstory canopy cover. STUDY AREAS Grand Portage Indian Reservation (Gichi Onigaming Ishkonigan) The Grand Portage Indian Reservation (Gichi Onigaming Ishkonigan in the Anishinaabe language) is located in the northeastern tip of Minnesota, USA (Figure 1). Gichi Onigaming encompasses approximately 22,600 ha of rugged terrain, ranging from 180 m at Lake Superior’s shoreline to 550 m on ridge tops. (Kraft et al. 2014). Climate and vegetation are significantly influenced by Lake Superior. Soils originated after the last ice age and con- sist of red rock outcrops, glacial moraines, and glacial lake superior clay plains (Kraft et al. 2014). Moose are present across the study area (0.27/km2), but their core range is inland away from Lake Superior (Oliveira-Santos et al. 2021). ALCES VOL. 60, 2024 TWIG PRODUCTION PATTERNS 79 Gichi Onigaming was about 80% for- ested and contained several forest cover types. In 2019, quaking aspen-paper birch (Betula papyrifera; Wiigwaas) forest repre- sented 56% of forest cover, with the remain- der composed of upland spruce-fir (Picea spp.- Abies balsamea; Mina’ig – Zhingob) forests (~12%), sugar maple (Acer saccha- rum; Ininaatig) forests (~6%), lowland swamp conifers (~5%), white cedar (Thuja occidentalis; Giizhik) swamps (~5%), and miscellaneous other forest types (Gichi Onigaming Forestry, unpublished 2019 GIS data). The Gichi Onigaming Forestry Department actively manage Reservation forests for multiple Tribal uses including commercial timber, firewood, moose habi- tat, and traditional gathering. Forests were managed using timber harvest, prescribed fire, tree planting, and a variety of tree tend- ing practices such as mechanical release, pre-commercial thinning, and pruning. Isle Royale National Park (Minong) Part of the ancestral homeland of the Gichi Onigaming Anishinaabeg is Minong (now known as Isle Royale National Park), an archipelago of islands in northwestern Lake Superior near the USA-Canada bor- der. The Gichi Onigaming Anishinaabeg are the historical occupants of Minong and have stewarded it for centuries. This rela- tionship was recognized in 2019 when the Federal Government designated Minong a “Traditional Cultural Property” of the Gichi Onigaming Anishinaabeg (U.S. National Fig. 1. Map of transect locations at Gichi Onigaming (Grand Portage Indian Reservation, Minnesota, USA) and Minong (Isle Royale, MI, USA) June–September 2019 with regional context and transect diagram (insets). TWIG PRODUCTION PATTERNS ALCES VOL. 60, 2024 80 Park Service 2019). Over the last decade, the Gichi Onigaming Anishinaabeg and Isle Royale National Park administration have established a working co-stewardship relationship. The climate and vegetation of Minong are heavily influenced by Lake Superior. Long ridges and valleys that run the length of the island also influence the vegetation and micro-climate. Forests are a mix of mesic sugar maple, yellow birch (Betula alleghaniensis; Wiinizik), and white cedar on the western end of the island, mixed paper-birch spruce-fir forests in the central part of the island, and spruce and fir boreal woodlands on the eastern end of the island and in a narrow band along Lake Superior (Sanders and Kirschbaum 2023). Moose have inhabited the island since the early 1900s and have influenced the vegetation in significant ways through browsing (De Jager et al. 2020). Moose density on Minong was higher than on the mainland during the study period. The humid influence of Lake Superior, Park Service regulations, and fire suppres- sion have limited the extent and severity of disturbances such as wildfire, timber har- vest, and development over the last 85 years. Although wildfires have occurred since 1940, they are not frequent or large (Kraft et al. 2010). Logging and development are prohibited by the Park Service. METHODS Sampling Design and Data Collection Model fitting data: Gichi Onigaming and Minong browse survey data. – We selected potential sampling site polygons in areas known to be used by moose within 400 m of a road or trail to increase survey crew effi- ciency, stratified to be evenly distributed across cover types, treated/untreated areas, and treatment types. At Minong, accessible moose use areas were stratified by ecosys- tem group (Nature Conservancy 1999). At Gichi Onigaming, accessible moose use areas were stratified into areas with recent management (“Treated”: timber harvest, shearing, planting release, or precommercial thinning) and those without (“Untreated”) using 2008–18 Gichi Onigaming treatment history maps (Gichi Onigaming Forestry, unpublished data). Untreated areas were then further stratified by cover type (decidu- ous forest, evergreen forest, mixed forest, shrub/scrub, woody wetlands) using US National Land Cover Database (Homer et al. 2015) and treated areas by treatment type. Once these strata were identified, at least 3 site polygons were chosen from each stratum using visual selection in ArcGIS Pro such that sites were well distributed across study areas, resembling systematic random sam- pling. Site polygons ranged from 0.8 ha to 319.7 ha (mean 36.3 ha) at Gichi Onigaming, and from 0.6 ha to 644.7 ha (mean 75.3 ha) at Minong. A minimum of 3 transect start points were visually placed within each site polygon. Start points were placed such that transects would not overlap. Field crews collected vegetation data at 297 transects across 97 sites (Table 1) between 6 June and 7 September 2019 in Gichi Onigaming, and between 15 July and 2 August 2019 on Minong. Of these, 150 tran- sects across 47 sites were at Gichi Onigaming, and 147 transects across 50 sites were on Minong (Fig. 1). At Gichi Onigaming, one very large treated site contained nine tran- sect start points, and three large treated sites received six transect start points each. Due to human error and difficult terrain, one Minong site had only two transects and another Minong site had only one transect. Map errors during stratification required site reclassification after data collection, result- ing in an unbalanced design (Table 1). ALCES VOL. 60, 2024 TWIG PRODUCTION PATTERNS 81 We established transects (20-m x 2-m) beginning at the start points, with azimuths determined by a random number (0-360) gen- erated on a cell phone application (Fig. 1). We recorded visual estimates of canopy cover to the nearest 5% at the center of each transect. Following Portinga and Moen (2015), plants on transects with 0–50% canopy closure were considered open canopy, and plants from tran- sects with 51–100% canopy closure were con- sidered closed canopy. For each tree within a transect, we recorded the tree species, total height, and total number of twigs. Due to time constraints, full measurements of shrubs were collected on only a subset of transects col- lected by one crew member. These transects were well distributed across possible sites and transects, resembling systematic random sampling. Of 3,191 observed shrubs, we mea- sured 542 shrubs (17%) across 127 transects on 62 sites (31 on Gichi Onigaming, 31 on Minong). We recorded shrub species, total height, and total number of twigs. Total num- ber of twigs for both trees and shrubs was defined as the count of all twig ends within the browse zone (0.5-2.5 m above the ground), including both intact twig ends and twig ends that had been recently browsed or stripped of leaves. Total height was visually estimated. Trees with foliage overhanging the transect plots were also recorded, even if no twigs were within the browse zone. Validation data: Gichi Onigaming thin- ning project monitoring data. – We used an independent, out-of-sample moose browse monitoring data set to assess the validity of our Table 1. Cover type and treatment status of moose browse study sites at Gichi Onigaming (Grand Portage Indian Reservation, Minnesota, USA) and Minong (Isle Royale, MI, USA) June–September 2019. Cover types and ecosystems reclassified per Severud et al. (2023) at Gichi Onigaming. Place Ecosystem Type Treated Cover Type n Gichi Onigaming Rock Outcrop No Rock Outcrop 2 Gichi Onigaming Shrub Wetland No Alder Swamp 3 Gichi Onigaming Boreal Forest No Mature Boreal Conifer 3 Gichi Onigaming Boreal Forest No Mature Boreal Mix 12 Gichi Onigaming Boreal Forest No Mature Aspen/Birch 2 Gichi Onigaming Maple Forest No Mature Maple Hardwood 3 Subtotal 25 Gichi Onigaming Boreal Forest Yes Boreal Clearcut 11 Gichi Onigaming Boreal Forest Yes Managed Conifer 10 Gichi Onigaming Maple Forest Yes Maple Clearcut 1 Subtotal 22 Minong Wetland Forest No White Cedar Swamp 2 Minong Boreal Forest No Mature Boreal Conifer 13 Minong Boreal Forest No Mature Boreal Mix 14 Minong Boreal Forest No Mature Aspen/Birch 7 Minong Maple Forest No Mature Maple Hardwood 14 Subtotal 50 Grand Total 97 TWIG PRODUCTION PATTERNS ALCES VOL. 60, 2024 82 twig production equations. These data were collected during fall 2021 to monitor moose browsing and vegetation response to an experi- mental pre-commercial aspen strip-thinning treatment conducted at Gichi Onigaming. Thinning sites were high productivity (aspen height at 50 yrs > 21.3 m [70 ft]) young aspen forests that had been clearcut 9-19 years prior and thinned in 2020 or early 2021. We collected browse monitoring data on 120 transects across 4 sites. At each site, we installed 10 transects in cut strips, 10 tran- sects in leave strips (uncut strips adjacent to cut strips), and 10 transects in control areas (contiguous areas without any treatment), for a total of 30 transects per site. Data were collected similarly to those described above, except that transects were 10 m x 2 m, treated area transects were oriented diagonally across strips, all shrubs were measured, and twigs counts over 15 were estimated to the nearest five twigs. Statistical Analysis Regression model forms. – We fit non-lin- ear regression models to our data to deter- mine which relationship with plant height and which species and canopy cover covari- ates best explained variation in accessible twig production, quantified by total twig count. Candidate model forms (Table 2) were a non-linear unimodal curve repre- senting Hypothesis 1 (unimodal model), an allometric growth curve representing the Grigal and Ohman equations (allometric model), and an intercept only model repre- senting the null hypothesis of no relation- ship between twigs and height (null model). We modeled the unimodal non-linear height relationship using a custom variant of the exponentiated exponential function (i.e. y = αλ(1-e-λx)α e-λx), which is used in sur- vival analysis as an alternative to Weibull and gamma distributions (Gupta and Kundu 2001, Aslam et al. 2010, Nadarajah 2011). For our model, we used an order α = 1 function, exponentiated the x-terms, and added normalizing and scaling parameters (Bolker 2008). Our parameterization took the form: y T e e4 s i s x H x H , ln2 2 ln2s i s c Bs Hs s i s c Bs Hs, ln 1 , ln 1 = −           ( ) ( )−       −       +       +       (1) where ys,i is the total twig count of observa- tion i of group s, Ts is a fitted coefficient for group s representing the mean peak (maxi- mum) twig count, xs,i is the height (in meters) of observation i of group s, Hs is a fitted coef- ficient for group s representing the mean height at which twigs peak, Bs is a fitted coefficient for group s representing the mean height increment above peak at which twig production declines by 90% (i.e. the recipro- cal of the rate at which twig production declines with height), c is a fixed normaliz- ing coefficient, in this case c = ln(5.284428) (Supplementary Material 1), ε(s,i) is the ran- dom error of observation i of group s, and s is a group determined by some combination of categorical covariates (e.g. species × can- opy). This parameterization allows biologi- cally meaningful comparisons between parameters of differing species or canopy status. It can be mathematically proven that the fitted parameters of equation 1 represent the location of the curve maxima at (H, T) and the height increment after the peak at which twig production declines by 90% (i.e. to 10% of maxima; H+B,0.1T) (Figure 2, Supplementary Material 1). For the allometric equation representing the Grigal and Ohman (1977) methodology, we used the equation: ε= +y A xs i s s i P s i, , , s (2) ALCES VOL. 60, 2024 TWIG PRODUCTION PATTERNS 83 where ys,i, xs,i , εs,i and s are as above, and As and Ps are a fitted allometric coefficients for group s. These coefficients are not com- parable to the coefficients of equation 1. Finally, the null model took the simple form ys,i = μs + εs,i where μs is the mean of group s and εs,i is the random error as above. We considered 3 potential error distributions (εs): Poisson, normal with standard deviation σs, and negative binomial with size parame- ter Ks. See below for error distribution selec- tion procedure. Regression model covariates. – To test hypotheses 2 and 3, we allowed model parameters related to magnitude (T, A, or µ) to vary according to species and canopy cover (open vs. shaded) covariates (group s) (Table 2). Options for magnitude parameter covariates were Species only, Canopy only, Species + Canopy, Species × Canopy, and intercept only. To prevent overfitting of spe- cies with small sample sizes, we only allowed curve shape and error distribution parameters to vary with either species or an overall intercept. We defined the fully satu- rated (i.e., all possible covariates) model as the model where magnitude parameters var- ied with the species x canopy interaction, and curve shape and error distribution parameters (if present) varied with species. Regression model fitting and selection. – All models were fit using the maximum likelihood mle2 function of the R package bbmle (Bolker et al. 2023), the package Table 2. Candidate models to estimate mean twig count within reach of moose per woody plant from height, species, and canopy cover at Gichi Onigaming (Grand Portage Indian Reservation, Minnesota, USA) and Minong (Isle Royale, MI, USA) June–September 2019. The B coefficient represents the mean height increment above peak to 90% decline in twig production within moose browse zone (i.e. the reciprocal of the rate of decline in twig production with height). Variables are defined as: ys,i is the total twig count of observation i of group s, xs,i is the height (in meters) of observation i of group s, c is a fixed normalizing coefficient, in this case c = ln(5.284428). Max(ys) is the maximum observed twig count of group s, and max(xs, 4) is the maximum observed height of group s or 4, whichever is greater. Possible random error distributions (εs) were Poisson, normal with standard deviation σs, or negative binomial with size parameter Ks. Model Equation Parameters Symbol Meaning Constraints Maximum Possible Covariates Unimodal y T e e 4 s i s x H x H , ln 2 2 ln 2 s i s c Bs Hs s i s c Bs Hs , ln 1 , ln 1 = −             ( ) ( ) −       −       +       +       (Eq.1) Ts mean peak twig count 0.1 ≤ T ≤ max(ys) Species x Canopy Hs mean height at peak count 0.5 m ≤ H ≤ max(xs) Species Bs mean height to 90% production decline 2 ≤ Bs ≤ max(4, xs) Species εs,i random error (Ks or σs) 0.001 ≤ σs ≤ max(ys) 0.01 ≤ Ks ≤ 100 Species Allometric ys,i = As * xs,i ^(Ps) + εs,i (Eq.2) As slope 0.001 ≤ As ≤ 100 Species x Canopy Ps power 0.001 ≤ Ps ≤ 20 Species εs,i random error (Ks or σs) 0.001 ≤ σs ≤ max(ys) 0.01 ≤ Ks ≤ 100 Species Null ys,i = μs + εs,i (Eq.3) μs mean twig count 0.001 ≤ μs ≤ max(ys) Species x Canopy εs,i random error (Ks or σs) 0.001 ≤ σs ≤ max(ys) 0.01 ≤ Ks ≤ 100 Species TWIG PRODUCTION PATTERNS ALCES VOL. 60, 2024 84 optimx (Nash 2014), and the L-BFGS-B bounded optimizing algorithm (Byrd et al. 1995). Models were evaluated using Akaike’s Information Criterion for small sample sizes (AICc; Burnham and Anderson 1998). Before fitting, we dropped species with fewer than 10 observations to prevent model non-convergence. We constrained model parameters between lower and upper bounds (Table 2) to ensure model conver- gence and keep parameters within observed values (Table 3). Notably, constraints of 2 and 4 for B prevented extreme values that crashed the statistical software. Fitting occurred in two stages following Zuur et al. (2009) and Bolker (2008): 1) error distribution fitting; 2) backwards selec- tion of covariates. In the first stage, we fit a fully saturated version of each of the 3 model forms with the 3 possible random error dis- tributions (Poisson, normal, or negative binomial). We then selected the best fit ran- dom error distribution for each of the 3 model forms using AICc. In the second stage, we performed backwards selection sepa- rately on each of the three fully saturated model forms with the error distribution selected in the first stage. Each possible Fig. 2. Interpretations of the parameters (T, H, B) and locations of critical points of Equation 1. Small B values indicate a high rate of production decline with height. Large B values indicate a low rate of production decline with height. ALCES VOL. 60, 2024 TWIG PRODUCTION PATTERNS 85 categorical covariate (or interaction) was sequentially removed from each parameter of each model form and the results evaluated by AICc. The variable removal with the greatest improvement in fit was selected, and the process repeated until no more vari- ables could be removed without worsening fit. If removing a variable did not worsen fit by more 2 AIC units, we considered that variable uninformative and removed it (Arnold 2010). Moose browsing preferences. – To aid interpretation of our results, we classified tree and shrub species by moose browse preferences reported by Portinga and Moen (2015; Tables 1, 8), realizing this is a coarse filter and many finer-scale variables influ- ence forage quality. For simplicity, we syn- thesized summer and winter preferences to classify overall preference. Species used sig- nificantly more often than available in at least one season were classified as preferred. Species not used significantly more often than available in at least one season and used according to abundance in at least one sea- son were classified as neutral. Species appearing in Portinga and Moen (2015; Tables 1, 8) and not classified as preferred or Table 3. Summary statistics of tree and shrub data used to fit models from moose browse study sites at Gichi Onigaming (Grand Portage Indian Reservation, Minnesota, USA) and Minong (Isle Royale, MI, USA) June–September 2019. Species Common Name Species Scientific Name n Height (m) Twigs Per Plant Med† σ Min Max Med† σ Min Max Balsam Fir Abies balsamea 2092 1.00 2.3 0.1 20 8 27.0 0 240 Red Maple Acer rubrum 24 1.90 1.7 0.5 5 24 20.1 4 102 Sugar Maple Acer saccharum 1217 0.60 1.9 0.1 18 4 10.1 0 110 Mountain Maple Acer spicatum 1709 0.80 1.1 0.1 8 6 11.0 0 114 Alder Spp. Alnus spp. 35 4.00 2.1 0.5 8 10 13.3 0 52 Juneberry Spp. Amelanchier spp. 83 1.20 0.7 0.3 3 5 11.2 0 46 Yellow Birch Betula alleghaniensis 197 1.30 2.8 0.5 15 8 10.1 0 91 Paper Birch Betula papyrifera 446 2.00 4.2 0.2 25 16 38.3 0 400 Beaked Hazelnut Corylus cornuta 66 1.00 0.9 0.2 6 12 18.4 0 122 Red Osier Dogwood Cornus sericea 59 1.20 0.2 0.4 1.2 5 2.5 1 21 Hawthorn Crataegus spp. 16 0.90 1.7 0.5 7 6 5.5 2 19 Bush Honeysuckle Diervilla lonicera 64 0.50 0.1 0.5 1 3 9.1 1 51 Black Ash Fraxinus nigra 249 1.50 3.3 0.2 18 14 19.4 0 142 Canada Honeysuckle Lonicera canadensis 27 1.00 0.4 0.5 1.8 6 13.0 3 70 White Pine Pinus strobus 20 2.00 2.2 0.2 8 30 47.9 0 150 Balsam Poplar Populus balsamifera 33 1.50 1.6 0.5 8 11 15.6 3 74 Quaking Aspen Populus tremuloides 948 2.00 3.7 0.1 20 8 15.3 0 113 Fire Cherry Prunus pensylvanica 62 1.25 0.9 0.3 5 13 20.8 0 100 Chokecherry Prunus virginiana 176 1.00 1.2 0.2 7 13.5 14.6 1 104 Willow Spp. Salix spp. 22 1.75 1.8 0.5 6 16 16.0 0 52 Red Elderberry Sambucus racemosa 10 1.50 0.9 0.6 4 12 5.1 1 18 Mountain Ash Sorbus spp. 776 0.70 0.7 0.1 7 4 5.0 1 47 White Cedar Thuja occidentalis 187 1.50 3.3 0.3 20 12 28.0 0 230 †Median TWIG PRODUCTION PATTERNS ALCES VOL. 60, 2024 86 neutral were classified as avoided/rarely used. Species not included in Portinga and Moen (2015; Tables 1, 8) were classified as not evaluated. Within sample model assessment. – We assessed in-sample accuracy and preci- sion by calculating the bias (mean error), RMSE (standard deviation of the error), Pearson correlation, and a pseudo-R2 (square of Pearson correlation) for the best fit model and leading alternative models. We calcu- lated approximate RMSE confidence inter- vals using the square root of variance confidence intervals [((n-1)*σ2)/χ2 (0.975,(n-1)) ≤ σ2 ≤ ((n-1)*σ2)/χ2 (0.025,(n-1))]. Differences in bias and RMSE between alternative models were assessed using repeated measures anal- ysis of variance (ANOVA), Mauchley’s test for sphericity, and post-hoc t-tests in R pack- age rstatix (Kassambara 2023). To control for the false discovery rate of multiple com- parisons, we used the YB correction (Benjamini and Yekutieli 2001). Production differences among species and canopy cover. – To understand twig production differences between species and canopy cover, we conducted post-hoc com- parisons on parameters from the final model using the glht and cldList functions in the multcomp and rcompanion R packages, respectively (Hothorn et al. 2023, Mangiafico 2024). To control for multiple comparisons, we used the fdr method (Benjamini and Hochberg 1995) for comparisons between species and the YB method for canopy com- parisons within species. Out-of-sample model validation. – To validate the model, we applied the final model equations to the independent Gichi Onigaming thinning project browse moni- toring data and evaluated out-of-sample model performance at the population level. Regression equations have at least 2 sources of error: the prediction error of individual observations and the estimation error of population parameters (Pardoe 2020). In applications such as monitoring concerned with population means, individual observa- tion errors are unimportant because oppos- ing errors cancel out. Therefore, we evaluated performance at the population level using summed transect data. We used open canopy equation coefficients for cut strips, closed canopy equation coefficients for leave strips and control areas. We evaluated the accuracy and precision of the best fit model and alter- native models using the same statistical measures as within-sample model assess- ment, as well as by fitting linear regression lines to the modeled vs. observed data. All data analyses were conducted in the R statistical environment version 4.3.3 (R Core Team 2024). We determined signifi- cance according to P < 0.05 and estimated 95% confidence intervals. Simulation of Forest Management Implications To understand the forest management impli- cations of our results, we coupled our moose forage equations with tree age-height curves to simulate forage dynamics over time for three tree species (quaking aspen, paper birch, and red maple (Acer rubrum; Zhiishiigimewanzh) known to be preferred by moose (Irwin 1985, Portinga and Moen 2015). For simplicity, we assumed pure even-aged stands with no understory shrubs and no tend- ing treatments. We used site index equations for the lake states (Michigan, Minnesota, Wisconsin) to model tree height growth as a function of age and site quality (Carmean et al. 1989). We assumed an aspen site index of 18.3 m (60 ft) at 50 years, typical of average sites in northeastern Minnesota, and calculated equiv- alent site indexes for other species using equa- tions in Carmean et al. (2013). We assumed initial stem densities and mean annual self-thinning mortality rates (Supplementary Material 2, Table S2.1) based on USFS ALCES VOL. 60, 2024 TWIG PRODUCTION PATTERNS 87 silviculture guides (Leak et al. 1969, Perala 1977, Safford 1983). We multiplied twig den- sities by estimates of average mass per twig (Portinga and Moen 2015) to calculate forage biomass densities for both winter and summer. We compared simulation results to a published minimum forage density threshold for moose use of 50 kg/ha (Allen et al. 1987). RESULTS We identified and modeled 23 species (11 trees and 12 shrubs) from 8,518 total obser- vations (Table 3). Trees and shrubs were identified to species-level except for willows (Salix spp., Oziisigobiminzh), juneberries (Amelanchier spp.; Gozigwaakominaga- awanzh; alders (Alnus spp.; Wadoop), haw- thorns (Crataegus spp.; Miinensagaawanzh), and mountain ash (Sorbus spp.; Makominaga- awanzh), which were identified to genus only. Best Fit Model: Relationship Between Twig Production and Covariates Negative binomial error distributions fit best across all 3 model forms (Table 4). The uni- modal model best described the relationship between twig production and woody plant height (Table 5). The top 3 unimodal models fit better than the best allometric or null models in all cases and by a large margin (Table 5). The saturated unimodal model fit best overall, indicating that peak number of twigs per plant (T) varied with the interaction between canopy and species, and all other parameters (H, B, K) varied by species (Figure 3, Table 5). The saturated model fit significantly better than the next closest model (Likelihood ratio test χ2 = 218.74, P < 0.001). The saturated unimodal model had a moderately strong correlation between observed and predicted data (0.58 ± 0.02), explained roughly a third of the observed variance, and had a higher correlation than the saturated allometric or null models (Supplementary Material 3, Table S3.1). RMSE varied significantly among models (Mauchley’s test w = 0.39, P < 0.001), and the unimodal model had the lowest RMSE among models (Table S3.1). Bias differed significantly among models (repeated measures ANOVA F = 10.5, adj.DFn = 1.93, adj.DFd = 16428.7, P < 0.001). Bias was not significant for the unimodal and null models (Supplementary Material 3, Table S3.2). The allometric model had sig- nificant over-estimation bias and signifi- cantly greater bias than the unimodal and null models. Unimodal model bias and RMSE varied among species (Table 6). Species-level bias ranged from + 0.99 for juneberry to -1.08 for red elderberry (Sambucus racemosa; Bibigwewanashk), but was not significant for any species. Table 4. Measures of fit (AICc) of 3 error distributions applied to 3 models used to estimate mean twig count within reach of moose per woody plant at Gichi Onigaming (Grand Portage Indian Reservation, Minnesota, USA) and Minong (Isle Royale, MI, USA) June–September 2019. All models fit with fully saturated parameter covariates (see text). Model Unimodal Allometric Null Error Distribution AICc AICc AICc Negative Binomial 55,287.7 56,725.6 59,594.7 Normal 66,396.9 68,186.2 69,853.1 Poisson 109,710.8 132,381.4 156,544.8 TWIG PRODUCTION PATTERNS ALCES VOL. 60, 2024 88 Of 23 modeled species, all coefficients were within bounds for 11 species, indicat- ing good predictive value (Table 6). Eight species had B coefficients (rate of decline with height) at an upper bound and other coefficients within bounds, indicating rea- sonable predictive value except at large heights (Table 6). The remaining 4 shrub species (red osier dogwood [Cornus sericea; Miskwaabiimizh], hawthorn, red elderberry [Sambucus racemosa; Bibigwewanashk], and bush honeysuckle [Diervilla lonicera]) had 2 coefficients at a bound, indicating lim- ited predictive or inferential value (Table 6). Of these, bush honeysuckle had a B coeffi- cient with an undefined standard error, indi- cating poor model convergence for this species. Differences in Peak Twig Production (T) Between Species and Canopy Cover Peak number of twigs was significantly greater in the open than under shade for 7 species (mountain maple, quaking aspen, mountain ash, balsam fir, black ash (Fraxinus nigra; Baapaagimaak), white cedar, yellow birch; Figure 4, Table 6). There was no sig- nificant difference in peak twig production between open and closed canopy for the other 16 species. Under open canopy conditions, mean peak number of twigs per plant (Topen) varied considerably among species, with seven spe- cies having greater peak production than aspen. In the open, 3 preferred species (paper birch, red maple, chokecherry (Prunus virgin- iana; Asasaweminagaawanzh) and 4 avoided Table 5. Measures of fit (AICc and likelihood ratio tests) for models and parameter covariates used to estimate mean twig count per plant within moose browse zone at Gichi Onigaming (Grand Portage Indian Reservation, Minnesota, USA) and Minong (Isle Royale, MI, USA) June–September 2019. Only the top three models for each equation form are displayed. K is the size parameter of the negative binomial distribution, df is number of model parameters. χ2, Δdf, and P-value are for the likelihood ratio tests for difference between top 3 models (1 vs. 2 and 1 vs.3). Equation Form Parameter Covariates Model Fit T H B K df AICc ΔAICc Chi-Sq Δdf P Unimodal Spp x Canopy Spp Spp Spp 115 55,288 0 Unimodal Spp + Canopy Spp Spp Spp 93 55,461 174 218.7 22 <0.001 Unimodal Spp x Canopy Spp 1 Spp 93 55,675 387 432.3 22 <0.001 A P K Allometric Spp x Canopy Spp Spp 92 56,726 1,438 Allometric Spp + Canopy Spp Spp 70 56,959 1,671 Allometric Spp x Canopy 1 Spp 70 57,022 1,735 µ K Null Spp x Canopy Spp 69 59,595 4,367 Null Spp + Canopy Spp 47 59,757 4,529 Null Spp x Canopy 1 47 60,382 5,153 ALCES VOL. 60, 2024 TWIG PRODUCTION PATTERNS 89 or not evaluated species (white pine, balsam fir, white cedar, black ash) had a significantly greater peak number of twigs than quaking aspen (Figure 4, Table 6). Of these, paper birch, white pine, and white cedar were the most productive. In the open, one preferred shrub (mountain ash) had a significantly lower peak number of twigs than quaking aspen. Under a shaded canopy, mean peak number of twigs per plant (Tclosed) also varied considerably among species. Six preferred species (paper birch, red maple, choke- cherry, fire cherry (Prunus pensylvanica; Bawa’iminagaawanzh), mountain maple, mountain ash), 3 avoided species (balsam fir, black ash, beaked hazel), and 2 not Fig. 3. Observed data and best fit model prediction of mean number of twigs per plant within moose browse zone as a function of height, species, and canopy cover at Gichi Onigaming (Grand Portage Indian Reservation, Minnesota, USA) and Minong (Isle Royale, MI, USA) June-September 2019. Point color intensity is proportional to number of observations. Canopy cover is open (0-50%) or shaded (51-100%). Equation is y = 4T(exp(-ln(2)(x/H)^(c/ln(1+B/H))) - exp(2ln(2)(x/H)^(c/ ln(1+B/H))) ), where y is mean number of twigs per plant, x is total height in meters, and coefficients are in Table 6. Symbols after species names indicate moose preference (+), proportional use (0), avoidance/rare use (-) or not evaluated (no symbol) in Portinga and Moen (2015). Observations beyond x = 20 and y = 120 not shown. TWIG PRODUCTION PATTERNS ALCES VOL. 60, 2024 90 Table 6. Coefficients (β) and standard errors (SE) and RMSE of best fit unimodal regression equation to predict mean number of twigs per plant within moose browse zone as a function of height, species, and canopy cover at Gichi Onigaming (Grand Portage Indian Reservation, Minnesota, USA) and Minong (Isle Royale, MI, USA) June–September 2019. Species T (Open) T (Shaded) ΔT H B K Error Common Name β SE β SE Δβ SE p β SE β SE β SE RMSE Bias Balsam Fir 48.3 2.3 40.2 2.2 8.1 2.0 0.003 4.7 0.1 12.5 0.6 1.1 0.0 23.1 -0.32 Red Maple 46.5 8.4 37.7 10.4 8.8 10.5 1.000 2.8 0.2 4.9 1.1 5.5 1.9 16.3 0.05 Sugar Maple 21.6 1.1 24.0 1.3 -2.4 1.1 0.206 3.2 0.1 7.0 0.3 2.6 0.1 8.0 -0.17 Mountain Maple 25.0 0.9 19.7 0.7 5.3 0.8 < 0.001 3.4 0.1 8.0# 0.5 3.7 0.2 8.6 -0.17 Speckled Alder 14.5 8.1 17.5 3.8 -2.9 9.4 1.000 4.5 0.9 8.0# 3.2 1.2 0.4 12.4 0.26 Juneberry Spp. 19.9 3.2 9.8 3.6 10.1 3.7 0.079 2.0 0.2 3.0# 0.5 1.4 0.3 11.1 0.99 Yellow Birch 24.4 3.7 14.8 1.7 9.6 2.9 0.016 3.0 0.2 6.5 0.6 2.5 0.3 8.9 0.17 Paper Birch 54.1 5.1 39.9 4.8 14.2 5.2 0.076 6.3 0.3 19.0 1.7 1.0 0.1 34.2 -0.98 Beaked Hazelnut 40.0 10.3 22.4 3.4 17.6 9.1 0.430 2.4 0.4 6.0# 1.5 2.0 0.4 17.3 0.51 Red Osier Dogwood 6.4 1.1 3.0 1.8 3.4 2.1 0.656 0.5# 0.2 3.0# 1.1 56.4 68.3 2.4 0.00 Hawthorn 12.6 2.6 19.0# 6.1 -6.4 6.2 1.000 2.5 0.3 7.0# 1.3 6.3 4.6 4.5 -0.78 Bush Honeysuckle 7.9 2.0 5.3 1.7 2.7 1.6 0.656 1.0# 0.3 2.0# NA 1.6 0.3 8.9 -0.09 Black Ash 44.9 4.3 31.5 2.7 13.4 3.7 0.009 4.3 0.2 10.8 0.8 2.4 0.2 16.9 0.47 Canada Honeysuckle 17.7 7.0 9.7 2.3 8.0 8.1 1.000 0.9 0.5 3.0# 1.1 1.7 0.5 13.3 0.74 White Pine 114.7 35.8 55.9 30.9 58.9 41.8 0.925 4.7 0.6 5.6 1.7 1.7 0.6 32.5 0.90 Balsam Poplar 31.1 7.9 24.6 7.9 6.4 11.1 1.000 3.9 0.6 8.0# 1.8 2.7 0.7 14.4 0.60 Quaking Aspen 23.9 1.1 9.9 1.2 14.0 1.5 < 0.001 2.4 0.1 4.5 0.2 1.1 0.1 13.0 -0.18 Fire Cherry 29.7 4.7 29.9 8.2 -0.2 8.3 1.000 2.5 0.4 5.0# 1.5 1.8 0.3 18.9 0.34 Chokecherry 36.6 3.6 33.9 2.7 2.8 3.1 1.000 3.1 0.2 6.2 0.8 4.3 0.6 11.5 0.16 Willow Spp. 44.3 16.6 29.9 9.8 14.5 19.2 1.000 2.4 0.3 3.2 0.5 2.4 0.9 13.2 0.65 Red Elderberry 18.0# 3.1 11.3 3.5 6.7 4.2 0.656 2.7 0.2 4.0# 1.2 66.0 77.7 2.4 -1.08 Mountain Ash 15.6 0.9 13.5 0.7 2.2 0.6 0.009 2.9 0.1 7.0# 0.6 8.8 1.1 3.7 0.00 White Cedar 66.3 11.0 35.1 4.4 31.2 9.6 0.018 3.2 0.1 5.3 0.3 1.4 0.2 23.5 0.66 Equation is y = 4T(exp(-ln(2)(x/H)^(c/ln(1+B/H))) - exp(-2ln(2)(x/H)^(c/ln(1+B/H)))) where y is mean number of twigs per plant, x is height in meters, T is the mean peak number of twigs per plant for the species and canopy combination, H is the mean height at which number of twigs peak, B is the mean height increment above peak to 90% decline in twig production (i.e. the reciprocal of the rate of decline in twig production with height), and c = ln(5.284428). K is the dispersion parameter of negative binomial error. Canopy cover is defined as open (0-50%) or shaded (51-100%). ΔT values are post-hoc comparisons corrected by Benjamini-Yekutieli (BY) method. Underlined ΔT values significant (p < 0.05). Underlined β values not significantly different from zero (p > 0.05). “#” indicates coefficient at boundary. “NA” indicates undefined. ALCES VOL. 60, 2024 TWIG PRODUCTION PATTERNS 91 evaluated species (white cedar, sugar maple) had significantly greater peak production than quaking aspen under canopy shade (Figure 4, Table 6). Species Differences in Height at Peak Production (H) Mean height at peak number of twigs per plant (H) varied considerably among species, with most species having a greater mean height at peak than aspen. Five preferred or used species (paper birch, mountain maple, chokecherry, mountain ash, willow), and 8 avoided or not evaluated species (white pine, balsam fir, speckled alder, black ash, balsam poplar, white cedar, sugar maple, yellow birch) had significantly greater mean height at peak than quaking aspen (Figure 5, Table 6). One preferred shrub (juneberry) and 1 not evaluated shrub species (Canada honey- suckle) had significantly lower height at peak than quaking aspen. Fig. 4. Modeled peak (T) of mean number of twigs per plant within moose browse zone by species and canopy cover at Gichi Onigaming (Grand Portage Indian Reservation, Minnesota, USA) and Minong (Isle Royale, MI, USA) June–September 2019. Symbols after species names indicate moose preference (+), proportional use (0), avoidance/rare use (-) or not evaluated (no symbol) in Portinga and Moen (2015). Letters indicate differences between species within the same canopy condition. Species not sharing any letter are different by Tukey-test with “fdr” correction (P < 0.05). Asterisks indicate differences between open and shaded conditions within species by Tukey-test with “YB” correction (*P < 0.05; **P < 0.01; ***P < 0.001). “#” indicates estimate at boundary. TWIG PRODUCTION PATTERNS ALCES VOL. 60, 2024 92 Species Differences in Rate of Twig Decline with Height (B) Among species with B coefficients within the range of observed values, one preferred tree (paper birch) and five avoided or not evaluated trees (balsam fir, black ash, sugar maple, and yellow birch) had significantly greater B coefficients than aspen, indicating lower rates of twig decline with height. One neutral shrub (willow) had a significantly lower B coefficient than aspen (Figure 6, Table 6), indicating a higher rate of twig decline with height. Among species with B coefficients at the boundary, indicating lim- ited inferential value, mountain maple and mountain ash had a greater B coefficient than aspen, whereas juneberry had a lower B coefficient than aspen. Out-of-Sample Model Validation at the Population Level The validation data set contained 18 of the 23 species fit in the original model, and bias and RMSE could be assessed for 17 species (Supplementary Material 4, Table S4.3). Unimodal estimates of out-of-sample number of twigs per transect explained 54.7% of the variance (Figure 7, Supplementary Material 4, Table S4.1). Correlations were higher for the unimodal model than the allometric and null models, although confidence intervals overlapped. Population level out-of-sample RMSE differed significantly among models (Mauchley’s test w = 0.249, p < 0.001). Unimodal model RMSE was greater than the null model and less than the allometric model. All models had significant out-of-sample overprediction bias at the population level (Figure 7, Supplementary Material 4, Table S4.2,). Bias differed significantly among models (repeated measures ANOVA F = 183.9, adj.DFn = 1.7, adj.DFd = 207.1, P < 0.001). The unimodal model had signifi- cantly higher bias than the null model and significantly lower bias than the allometric model. Out-of-sample bias and RMSE varied between species (Table S4.3). Overestimation bias was significant for 12 of 17 species. Simulation Model Results Simulation results for pure even-aged stands of 3 preferred forage species on average sites indicated that whereas quaking aspen pro- duced ~1.4-2.2 times more peak forage bio- mass than paper birch or red maple, paper birch and red maple produced usable forage densities for ~1.2-2.9 times longer (Table 7, Figure 8). In summer, quaking aspen had the highest, red maple the middle, and paper birch the lowest peak forage density. Summer paper birch forage was above the usable den- sity threshold for 26 years, red maple for 15 years, and aspen for 10 years. In winter, quak- ing aspen had the highest, paper birch the middle, and red maple the lowest peak forage density. Winter paper birch forage was above the usable density threshold for 26 years, red maple for 11 years, and aspen for 9 years. DISCUSSION This study and its implications support the cultural and ecological stewardship prior- ities of the Gichi Onigaming Anishinaabeg to advance moose restoration. This study’s spe- cies-specific twig production models offer a tool to evaluate and anticipate how different forest treatments influence moose forage dynamics over time, thereby supporting Indigenous-led strategies that seek to balance biodiversity, habitat health, and cultural con- tinuity. Moose are a key species ecologically and foundation of Indigenous food sover- eignty, ceremonial practice, and intergenerational knowledge transmission (Stults et. al 2016, Garwood et al. 2023, Moore et al. 2024, Severud et al. 2023). Forest management decisions involving timber har- vest and post-disturbance regeneration have significant implications for the viability of culturally essential moose populations. ALCES VOL. 60, 2024 TWIG PRODUCTION PATTERNS 93 Differences in Available Twig Production Between Species We have shown that production of twigs available to moose varies considerably among tree and shrub species (H2); that twig produc- tion is sometimes greater under open cano- pies (H3); that twig production increases, peaks, then declines as stem height increases (H1); and that the rate of decline with height varies somewhat among species (H1 and H2). Paper birch and balsam fir generally had the greatest peak twig production, greatest heights at peak production, and slowest rates of production decline. In contrast, quaking aspen generally had lower peak twig produc- tion, lower height at peak production, and faster rates of decline with height than most other trees and some tall shrubs. Two large shrubs (mountain maple, fire cherry) often preferred by moose (Miquelle and Jordan 1979, Belovsky 1981, Irwin 1985, Ricard and Doucet 1999, Portinga and Moen 2015) also performed favorably compared to quaking aspen. Mountain maple had similar peak twig production under open canopies, greater peak production under closed canopies, and greater height at peak production compared to aspen. Fire cherry had similar peak twig production Fig. 5. Modelled mean height (H) at peak number of twigs per plant within moose browse zone by species at Gichi Onigaming (Grand Portage Indian Reservation, Minnesota, USA) and Minong (Isle Royale, MI, USA) June-September 2019. Species not sharing any letter are significantly different by Tukey-test (P < 0.05) with “fdr” correction. “#” indicates estimate at boundary. Symbols after species names indicate moose preference (+), proportional use (0), avoidance/rare use (-) or not evaluated (no symbol) in Portinga and Moen (2015). TWIG PRODUCTION PATTERNS ALCES VOL. 60, 2024 94 under open canopies, greater peak production under closed canopies, and similar heights at peak production compared to aspen. This suggests that these shrubs, with typically slower height growth than aspen and ability to persist in the understory, have the potential to be significant and long persisting sources of preferred moose browse. Our twig production model appears rea- sonably consistent with observed growth- form characteristics of many tree and shrub species. The low peak twig production and rapid decline in twig production of quaking aspen among deciduous trees is consistent with findings that aspen have strong apical dominance (Farmer 1962), tend toward ver- tical rather than horizontal canopy exten- sion when young (Harper 2008) and self-prune lower branches more rapidly than other species (Puettmann and Reich 1995). In contrast, conifers such as balsam fir and white cedar have conical crown structure Fig. 6. Height to 90% production decline (B) by species at Gichi Onigaming (Grand Portage Indian Reservation, Minnesota, USA) and Minong (Isle Royale, MI, USA) June-September 2019. The B coefficient represents the mean height increment above peak to 90% decline in twig production (i.e. the reciprocal of the rate of decline in twig production with height). Species not sharing any letter are significantly different by Tukey-test (p < 0.05) with “fdr” correction. “#” indicates estimate at boundary. Bush honeysuckle excluded because estimate had undefined standard error. Symbols after species names indicate moose preference (+), proportional use (0), avoidance/rare use (-) or not evaluated (no symbol) in Portinga and Moen (2015). ALCES VOL. 60, 2024 TWIG PRODUCTION PATTERNS 95 and tend to retain their lower branches. White pine also has a conical crown and tends to delay self-pruning until at least 15 years old, after which it may self-prune aggressively (Burns and Honkala 1990). Moreover, evergreen conifer needles live longer than deciduous leaves. Theoretically, these factors together should result in more twigs within the browse zone. The higher peak twig production for these species in our model is consistent with this observa- tion (Figure 3, Table 6). Similarly, the slow decline in twig production of balsam fir rel- ative to other trees in our model (Figure 5) is consistent with the field observation that, of these three conifers, balsam fir tends to retain its lower branches the longest, retain- ing some lower twigs even at very tall heights. This may be in part due to the extremely high shade tolerance of balsam fir. Other studies have found that shade tol- erant evergreens tend to maximize leaf lon- gevity to survive under low light conditions (Walters and Reich 1999). Although somewhat surprising for an early successional species, the high peak twig production and slower decline in twig production of paper birch is also consistent with the known biology of paper birch. When grown on the same site, paper birch generally have slower height growth than quaking aspen (Carmean et al. 2013), sug- gesting less allocation of resources to height growth. Moreover, paper birch also produce abundant stump and epicormic sprouts when stressed or harvested (Burns and Honkala 1990). Based on field observations, we fur- ther speculate that although paper birch self- prune (Zenner and Puettmann 2008), they may do so more slowly than aspen. Fig. 7. Comparison of predictions by the 3 alternative models with out-of-sample (validation) observations of number of twigs per transect within moose browse zone at Gichi Onigaming (Grand Portage Indian Reservation, Minnesota, USA), September 2021. Black line is 1:1 line. Gray line is the predicted regression line, with regression equation and coefficient of determination listed on plot. TWIG PRODUCTION PATTERNS ALCES VOL. 60, 2024 96 The production of twigs available to moose may also be augmented by the effects of browsing itself. Previous research at Minong and in Sweden has shown that birch and aspen respond to light to moderate moose browsing by increasing the production of fine branches and foliage (De Jager and Pastor 2008, De Jager et al. 2009, Pastor 2016). For deciduous species with slower height growth, such as paper birch or red maple, this effect might amplify available twig production over time relative to faster growing species. Consistent with this explanation, our study observed greater peak twig production for both paper birch and red maple (Figure 1, Table 6). Interestingly, birch also had greater variability in twig production than aspen, despite a smaller sample size (Table 3, Figure 1). Although other explanatory mechanisms are possible, such a pattern might result if birch responds more strongly than aspen to variabil- ity in browsing pressure across the landscape. Differences in Twig Production with Canopy Cover Our finding that trees and shrubs generally produced more twigs in the open than under shade and that the difference varies between species is generally consistent with a large body of tree physiology studies and prior moose forage research (see reviews by Lusk et al. 2008 and Walters and Reich 1999). In general, the tree physiology literature finds that shade grown woody plants tend to reduce photosynthetic capacity (e.g. number and size of leaves) relative to open grown plants of the same species in order to man- age the balance between photosynthesis and respiration (Lusk et al. 2008). Further, dif- ferences in allocation in full sun vs. shade between species are related to shade toler- ance, evergreen vs. deciduous leaves, leaf longevity, and leaf mass (Walters and Reich 1999). An earlier study of moose browse found that mean moose bite mass is often – but not always – larger for open-grown plants than for shaded plants (Portinga and Moen 2015). Although size may also be related to other factors such as moose selec- tivity or optimal foraging, the Portinga and Moen (2015) result is consistent with the idea that shaded plants allocate less mass to leaves and twigs. Taken together, these bod- ies of work suggest that our model is biolog- ically reasonable. Model Validation and Limitations We have also shown that unimodal twig pro- duction-height equations perform better than allometric equations, and that our equations performed reasonably well with independent out-of-sample data as quantified with some of our validation metrics. Unimodal equa- tions consistently had higher pseudo-R2, lower bias, and lower RMSE than allometric equations on both fitting and validation data. The unimodal model was also consistently better correlated with actual twig counts Table 7. Comparison of simulated forage density (kg/ha) from ages 0-30 years between paper birch, quaking aspen, and red maple in winter and summer for aspen site index of 18.3 m (60 ft) at 50 years in northeastern Minnesota using best-fit unimodal model. Species Summer Winter First Year ≥ 50 kg/ha Last Year ≥ 50 kg/ha Max. Density (kg/ha) First Year ≥ 50 kg/ha Last Year ≥ 50 kg/ha Max. Density (kg/ha) Paper Birch 2 27 131 2 27 131 Quaking Aspen 1 10 294 1 9 189 Red Maple 2 16 195 3 13 98 ALCES VOL. 60, 2024 TWIG PRODUCTION PATTERNS 97 than the species-cover average model for both fitting and validation data. These results demonstrate the greater predictive ability of unimodal models over Grigal and Ohmann (1977) style allometric models. However, although our unimodal model did a good job reflecting differences in forage density between plots, it appeared to have a signifi- cant bias that consistently overestimated the absolute forage density of the validation data. This suggests that although the model was adequate for monitoring relative differ- ences between sites or changes in sites over time, it should be used with caution in appli- cations requiring high accuracy such as cal- culating moose carrying capacity. The model’s apparent lack of accuracy and ability to explain only 54.7% of variabil- ity is likely related to limited sample sizes, poor sampling balance, protocol deficien- cies, missing covariates and biological pro- cess variance. Several shrubs and trees had less than 50 samples, which limited the pre- cision of parameter estimates. For these spe- cies, fitting data was poorly balanced with respect to height and site factors such as can- opy cover, productivity, and past manage- ment. For trees, small sample sizes were due to species rarity. For most shrubs, small sam- ple sizes were due to our subsampling of twig counts. Our measurement of multi- stemmed plants as one individual likely increased count variability and reduced esti- mate precision. Lack of identification to spe- cies for some shrubs may have pooled counts across species with divergent growth forms. For example, Salix discolor reaches heights of 10 m, whereas Salix humilis only reaches Fig. 8. Comparison of simulated forage density (kg/ha) from ages 0-30 years between paper birch, quaking aspen, and red maple in winter and summer for aspen site index of 18.3 m (60 ft) in northeastern Minnesota using best-fit unimodal model. Solid gray line indicates winter density. Dashed black line indicates summer density. Dotted line represents minimum forage density for moose use of 50 kg/ha. Paper birch winter and summer curves overlap because biomass/twig is the same. TWIG PRODUCTION PATTERNS ALCES VOL. 60, 2024 98 heights of 3 m (Smith 2008). Other covari- ates that might also influence twig produc- tion, such as stem diameter, neighboring stem density, browsing pressure, and abiotic factors affecting site productivity were not explicitly included in the model. Conversely, apparent bias might have resulted from an unrepresentative set of validation data; vali- dation data were collected at sites that all had above average site productivity (aspen height at 50 yrs > 21.3 m (70 ft)) and high stem densities resulting from past clearcut- ting. Compared with unmanaged sites with average productivity, the validation sites might have had above average rates of branch senescence from rapid height growth and heavy competition. It is possible our model may perform better or underestimate twig production elsewhere. Model accuracy could likely be improved by collecting additional fitting data with improved methods. In the future, twigs should be counted on all shrubs, rather than subsampling. Sampling should inten- tionally seek out species used by moose that were under-sampled in this study and better identify shrubs (e.g., willow, juneberry) to species or species groups. Improved proto- cols for dealing with multi-stemmed plants should be developed. Accounting for addi- tional factors such as neighboring stem den- sity and site productivity should be explored. Using Models with Forestry Field Surveys to Estimate Moose Forage Density Despite the limitations of our unimodal model noted above, the high correlation with validation data and consistency with pub- lished literature suggest the model is a useful first approximation to estimate twigs avail- able to moose for forestry and wildlife mon- itoring purposes. Applying these equations to common forest regeneration survey data (e.g. plot tallies of tree and shrub stems by species and height class) and multiplying by published twig mass values, a practitioner can estimate moose forage densities (e.g., kg per hectare or lbs per acre) by species. Forage densities can then be further grouped by moose preference (and avoidance) to esti- mate total acceptable browse. We developed a detailed protocol and worked example of the application of these equations to forest regeneration survey data (Supplementary Material 5, Table S5.1). The tally data used in that exercise were slightly modified from an actual plot taken during a routine regeneration survey in a dense part of a 6-year-old clearcut at Gichi Onigaming. Our estimate of 595 kg/ha, although not rep- resentative of the overall stand due to very high stem density on this plot, was similar to other high forage density estimates in the lit- erature (Crête 1989; Peterson et al. 2022). Our experience suggests that estimating for- age availability by applying our equations to already occurring forest regeneration sur- veys may be more efficient than conven- tional forage transect methods. It took less than 14 minutes to collect the forest regener- ation plot used in this example (Table S5). In contrast, it took on average ~60 minutes to collect the 2 m x 20 m forage transects used for fitting and ~30 minutes to collect the 2 m x 10 m transects used for validation. This suggests a more efficient use of field time, albeit with lower accuracy. Implications for Moose and Forest Management If abundance of diverse, high-quality forage can influence moose at the population level (Schrempp et al. 2019, Felton et al. 2020, Hayes et al. 2022, Peterson et al. 2022), then our results indicate that the assumption that aspen management alone is sufficient to main- tain moose populations should be re-exam- ined. To the extent that moose preference is a proxy for nutritional quality, our finding that ALCES VOL. 60, 2024 TWIG PRODUCTION PATTERNS 99 quaking aspen’s low twig production per stem and rapid rate of decline with height relative to other preferred tree and shrub species suggests that aspen alone may not provide a sufficient forage base for moose over time. Our simula- tion further supports the idea that quaking aspen’s rapid decline in twig production may scale up to the stand level, with consequences for moose forage availability on the landscape over time. Although aspen in our simulation produced a huge pulse of dense forage, the pulse was short-lived (9-10 years); paper birch and red maple provided forage for longer peri- ods of time (11-25 years). This illustrates an important paradox: although quaking aspen is preferred by moose (Irwin 1985, Portinga and Moen 2015), commercially valuable, and easy to regenerate by clearcutting, each stem does not provide many twigs, nor do they provide twigs once they have grown beyond a certain height. Although aspen’s low twig production per stem is offset by high stem densities in clearcuts, aspen’s rapid early height growth and self-pruning soon moves twigs out of the reach of moose. Conversely, our simulation suggested that young red maple and paper birch, two species preferred by moose but with low commercial value in northeastern Minnesota, may provide moose forage for longer periods of time than aspen. Managers in northeastern Minnesota have not historically managed specifically for paper birch due to low stumpage values and regeneration difficulties. More often, paper birch has been viewed as an incidental com- panion species and sometimes as a weed com- peting with more commercially valuable species. Red maple, less common and not typ- ically found in pure stands in northeastern Minnesota, has similarly often been viewed as a competing species to be eliminated. However, our results suggest that maintaining these tree species on the landscape, despite their low commercial value, may be important to sus- taining moose populations. Similarly, our results indicating favor- able performance of fire cherry and mountain maple relative to aspen suggest that maintain- ing preferred shrubs on the landscape may also be important to sustaining moose popula- tions. Foresters have historically viewed shrubs as a barrier to tree establishment, using herbicides and mechanical release to reduce shrub density. Our results suggest that a shift in perspective from “shrubs as competition” to “preferred shrubs as wildlife asset” may benefit foresters managing for moose. Such a shift would be consistent with the Anishinaabeg understanding that all species have value (Davidson-Hunt et al. 2005). Similarly, our results suggest that stands with a mix of preferred deciduous tree and shrub species may benefit moose more than stands with low species diversity. If, as our simula- tion found, different species produce peak forage at different times, then mixed stands would have several forage peaks resulting in more sustained forage availability. This might in turn benefit moose populations through improved calving success, associated with high preferred browse density (Hayes et al. 2022), and/or enhanced calf weight gain, which has been associated with higher decid- uous browse diversity (Felton et al. 2020). Although our results also suggest that balsam fir may provide abundant moose for- age for long periods of time, we caution that moose appear to avoid eating fir when other options are available (McNicol and Gilbert 1980, Cumming 1987, Risenhoover 1987, Hodgson 2010, Portinga and Moen 2015). Although high density moose populations on Minong often browse fir heavily in winter (Brandner et al. 1990), studies on Minong found that fir is eaten less than available (Hodgson 2010), browsed at lower rates than deciduous species (Risenhoover 1987), or are browsed more than cedar but less than decid- uous species (Belovsky 1981). Similarly, sev- eral studies of moose in Minnesota and TWIG PRODUCTION PATTERNS ALCES VOL. 60, 2024 100 northwest Ontario found that moose eat fir less than available (McNicol and Gilbert 1980, Cumming 1987, Portinga and Moen 2015), although one study in central Ontario found moose preferred fir (Routledge and Roese 2004). Possible reasons for avoidance are that fir has lower digestibility than decid- uous browse (Risenhoover 1987, Parikh et al. 2016) and contains terpenes (Terra-Berns 1993) that may be energetically costly for moose to detoxify (Parikh et al. 2016; but see also Hoy et al. 2022). Overall, our results suggest the impor- tance of managing intentionally for a diver- sity of preferred deciduous tree and shrub species in areas intended for moose forage, such as timber harvests. Post-harvest treat- ments such as mechanical site preparation, burning, shearing, and mechanical release could be used to increase the abundance, diversity, and longevity of preferred moose forage species. Mechanically releasing paper birch from aspen can shift stand composition towards birch (Zenner and Puettmann 2008). Light shearing or brush-sawing can also stim- ulate vigorous re-growth of red maple, moun- tain maple, fire cherry, aspen, and paper birch from stump or root sprouts (Krefting et al. 1956, Scheiner et al. 1988, Burns and Honkala 1990). Fire cherry and paper birch regenerate well from seed following more severe mechanical site preparation or prescribed fire (Burns and Honkala 1990). At least one Anishinaabeg community in what is now Canada historically burned specifically to improve moose forage (Theriault 2006), and many Tribal nations elsewhere in North American burn to increase forage for ungu- lates (Kimmerer and Lake 2001, Aguilar 2005, Connor et al. 2022, Hoagland and Albert 2023). Using prescribed fire to improve moose forage would support restoration of Anishinaabeg cultural burning, which was widely practiced historically (Larson et al. 2021, Johnson et al 2022). Conversely, herbicide use for release or site preparation should be avoided were possible because it can reduce the density of preferred deciduous browse species (Lautenschlager 1992, Johnson and Rea 2024). Other methods should be evaluated for impact on preferred forage density over time. Although promising, managing for pro- longed availability of preferred moose for- age rests on the assumption that moose prefer forage of high nutritional quality. Although other studies have drawn links between forage quality and moose prefer- ence (e.g. Belovsky 1981, Peterson et al. 2020), our study does not address either preference or quality. Measurement of for- age quality and moose preference is also notoriously complex. Forage quality (e.g. energy, protein, and mineral content) varies considerably across plant parts, species, sea- sons, and geographic regions, and is subject to complicating factors such as tannins and secondary metabolites that can limit uptake (Robbins et al. 1987, Lindroth 2001, Wam et al. 2018). Methods in moose preference studies vary considerably, often lacking quantitative rigor (e.g., Peek et al. 1976) or using measures such as Ivlev’s electivity (e.g. Irwin 1985) which may be biased and unreliable (Manly et al. 2002). As noted in the case of balsam fir, conflicting results between preference studies are common. More research is needed to better understand moose forage quality and its relationship to moose forage preferences in the western Lake Superior region. Managing for prolonged forage avail- ability may also have downsides. Quaking aspen can sometimes be browsed so heavily and repeatedly that it cannot mature, to the detriment of species reliant on mature and old aspen (Mathisen et al. 2017, Painter et al. 2018). Prolonged use in concentrated areas might also facilitate parasite transmission (Ellingwood et al. 2020, Hoy et al. 2021). ALCES VOL. 60, 2024 TWIG PRODUCTION PATTERNS 101 Care must be taken to ensure an abundance of forage across the landscape to prevent animal concentration and over browsing. Post-harvest treatments such as burning, shearing, and/or mechanical release may also reduce parasites such as brainworm (Parelaphostrongylus tenuis) and winter tick (Dermacentor albipictus) implicated in the decline of moose in northeastern Minnesota (Wünschmann et al. 2015, Carstensen et al. 2018, Severud et al. 2019, 2022; Oliveira- Santos et al. 2021). Prescribed fire can reduce densities of some slug and snail (gas- tropod) species that transmit brainworm to moose, while increasing densities of others (Nekola 2002). Prescribed burning can also reduce winter tick densities when litter and duff layers are sufficiently consumed (Drew et al. 1985), although winter tick populations might also rebound if enhanced forage increases moose abundance and density (Ellingwood et al. 2020, Hoy et al. 2021). At Gichi Onigaming, forestry treatments within the last 5 years (including mechanical release only) reduced overall gastropod densities and specifically densities of Deroceras spp., a frequent brainworm host (Severud et al. 2023). Finally, our results indicate that canopy openness is critically important to stimulat- ing high density moose forage. Managers should be mindful of light conditions as we seek to improve moose forage. Reserve areas in timber harvests, although import- ant for winter snow refuge (Mastenbrook and Cumming 1989) and summer cooling (Street et al. 2016), should not be so close together that they create excessive shade. In some mixed wood stands, early, heavy, and frequent thinning might also be used to maintain high densities of preferred forage over time while retaining mature trees for cover. In summary, our findings suggest that young, open canopied forests with a diverse mix of deciduous trees and shrubs likely provide abundant moose forage for longer durations than do the nearly pure aspen for- ests that often follow clearcutting alone in northeastern Minnesota. Species used and preferred by moose such as paper birch, red maple, mountain maple, and fire cherry may be of particular importance due to equal or greater peak forage production and/or slower rates of forage decline. In many aspen-dominated forests, promoting tree and shrub diversity to sustain adequate moose forage will likely require additional post-har- vest treatments such prescribed fire, site preparation, winter shearing, or mechanical (brush saw) release. These treatments may also benefit moose by reducing parasite den- sities on the landscape. Our findings reinforce the value of het- erogeneous post-disturbance regeneration in supporting sustained moose forage avail- ability, especially when guided by spe- cies-specific forage production curves. These outcomes align closely with Indigenous co-stewardship principles that prioritize ecological resilience, interspecies relationships, and culturally grounded resource governance (Moore et al. 2024). Collaborative research such as this study exemplifies how Indigenous science and western methodologies can be integrated to inform adaptive management. Future appli- cations of these equations, when paired with ongoing monitoring of moose health, preda- tion, and parasitism pressures (Wolf et al. 2021, Garwood et al. 2023, Severud et al. 2023, Weesies, in review) serve as a power- ful co-management tool. They offer the potential to tailor silvicultural prescriptions in a way that supports both moose recovery and food sovereignty goals of the Gichi Onigaming Anishinaabeg and affirm the vital role of Indigenous leadership in shap- ing forest futures that are ecologically func- tional and culturally enduring. TWIG PRODUCTION PATTERNS ALCES VOL. 60, 2024 102 ACKNOWLEDGEMENTS Funding for this project was provided by the United States Fish and Wildlife Service Tribal Wildlife Grant F19AP00035. We thank the Gichi Onigaming Tribal Council for their ongoing commitment to moose stewardship and research, and supervisors T. Miller and V. Cook for supporting this work. Many thanks to K. Woerheide, E. J. Isaac, T. Garwood, T. Walters, S. Young, and K. Gallup for collecting data. H. Fox provided invaluable GIS and data management sup- port. We thank E. Redix (Lac Courte Oreilles Band of Lake Superior Chippewa) for review and assistance with correct usage of Anishinaabeg language. We thank forestry technician and tribal moose hunter R. Spry for generating research questions and shar- ing his understanding of moose habitat use, and E. Carlson for sharing his understanding of Anishinaabeg fire use. We thank the Intertribal Timber Council for help in under- standing indigenous fire use in North America more broadly. N. DeCesare and 2 anonymous reviewers provided helpful feed- back that greatly improved this manuscript. This research leveraged data gained from a long-term ecosystem health research pro- gram led by the Gichi Onigaming Anishinaabeg (Grand Portage Band of Lake Superior Chippewa) and University of Minnesota. SUPPLEMENTARY MATERIAL Supplementary data are available with the article at https://alcesjournal.or/indexphp/ alces/article/view/1961. LITERATURE CITED AguilAr, george W. Sr. 2005. When the river ran wild! Indian traditions on the Mid-Columbia and the Warm Springs Reservation. University of Washington Press, Seattle, USA. Allen, A.W., P.A. JordAn and J.W. Terrell. 1987. Habitat suitability index models: moose, Lake Superior region. Biological Report 82(10.155). U.S. Department of the Interior, Fish and Wildlife Service, Washington, D.C. USA. Arnold, T.W. 2010. Uninformative parame- ters and model selection using Akaike’s information criterion. The Journal of Wildlife Management 74:1175-1178. AslAm, m., d. Kundu and m. AhmAd. 2010. Time truncated acceptance sampling plans for generalized exponential distri- bution. Journal of Applied Statistics 37:555-566. BelovsKy, g.e. 1981. Food plant selection by a generalist herbivore: The moose. Ecology 62:1020-1030. BenJAmini, y. and y. hochBerg. 1995. Controlling the false discovery rate: a practical and powerful approach to multi- ple testing. Journal of the Royal Statistical Society. Series B (Methodological) 57: 289-300. BenJAmini, y. and d. yeKuTieli. 2001. The control of the false discovery rate in multiple testing under dependency. The Annals of Statistics 29:1165-1188. BolKer, B.m. 2008. Ecological models and data in R. Princeton University Press, Princeton, N.J. BolKer, B., r.d.c. TeAm and i. giné- vázquez. 2023. bbmle: Tools for General Maximum Likelihood Estimation. (ver- sion 1.0.25.1). https://cran.r-project.org/ web/packages/bbmle/index.html. BrAndner, T.A., r.o. PeTerson and K.l. risenhoover. 1990. Balsam fir on Isle Royale: Effects of moose herbivory and population density. Ecology 71:155-164. BurnhAm, K.P. and d.r. Anderson. 1998. Practical Use of the Information-Theoretic Approach. Pages 75-117 in K. P. Burnham and D. R. Anderson (editors). Model Selection and Inference: A Practical Information-Theoretic Approach. Springer, New York, NY. USA. https://alcesjournal.or/indexphp/alces/article/view/1961 https://alcesjournal.or/indexphp/alces/article/view/1961 https://cran.r-project.org/web/packages/bbmle/index.html https://cran.r-project.org/web/packages/bbmle/index.html ALCES VOL. 60, 2024 TWIG PRODUCTION PATTERNS 103 Burns, r.m. and B.h. honKAlA. 1990. Silvics of North America. U.S. Department of Agriculture Forest Service, Agriculture Handbook 654. Byrd, r.h., P. lu, J. nocedAl and c. zhu. 1995. A Limited Memory Algorithm for Bound Constrained Optimization. SIAM Journal on Scientific Computing 16:1190-1208. cArmeAn, W.h., J.T. hAhn and r.d. JAcoBs. 1989. Site index curves for forest tree species in the eastern United States. U.S. Department of Agriculture Forest Service, Northern Research Station, General Technical Report NC-128. cArmeAn, W.h., J.T. hAhn, r.e. mcroBerTs and d. KAisershoT. 2013. Site index comparisons for forest species in the Upper Great Lakes area of the United States and Canada. U.S. Department of Agriculture Forest Service, Northern Research Station, General Technical Report NRS-113. cArsTensen, m., e.c. hildeBrAnd, d. PlATTner, m. dexTer, A. WünschmAnn and A. Armien. 2018. Causes of non-hunting mortality of adult moose in Minnesota, 2013-2017. Pages 30–40 in L. Cornicelli, M. Carstensen, M. Larsen, N. Davros, and B. Davis, editors. Summaries of wildlife research findings, 2017. Minnesota Department of Natural Resources, Wildlife Populations and Research Unit, St. Paul, Minnesota, USA. crêTe, M. 1989. Approximation of K carry- ing capacity for moose in eastern Quebec. Canadian Journal of Zoology 67:373-380. connor, T., e. TriPP, B. TriPP, B.J. sAxon, J. cAmArenA, A. donAhue, d. sArnA- WoJcicKi, l. mAcAulAy, T. BeAn, A. hAnBury-BroWn, and J. BrAshAres. 2022. Karuk ecological fire manage- ment practices promote elk habitat in northern California. Journal of Applied Ecology 59:1874-1883. cumming, H.G. 1987. Sixteen years of moose browse surveys in Ontario. Alces 23:125-156. dAvid, A.J., J.c. zAsAdA, d.W. gilmore and s.m. lAndhäusser. 2001. Current trends in the management of aspen and mixed aspen forests for sustainable production. The Forestry Chronicle 77: 525-532. dAvidson-hunT, i.J., P. JAcK, e. mAndAmin and B. WAPioKe. 2005. Iskatewizaagegan (Shoal Lake) Plant Knowledge: An Anishinaabe (Ojibway) Ethnobotany of Northwestern Ontario. Journal of Ethnobiology 25:189-227. de JAger, n.r. 2008. Multiple scale spatial dynamics of the moose-forest-soil eco- system of Isle Royale National Park, MI, USA. Ph.D. Thesis. University of Minnesota. 148 pp. de JAger, n.r. and J. PAsTor. 2008. Effects of moose Alces alces population density and site productivity on the canopy geometries of birch Betula pubescens and B. pendula and Scots pine Pinus syl- vestris. Wildlife Biology 14:251-262. de JAger, n.r., J. PAsTor and A.l. hodgson. 2009. Scaling the effects of moose browsing on forage distribution, from the geometry of plant canopies to land- scapes. Ecological Monographs 79: 281-297. de JAger, n.r., J.J. rohWeder and m.J. duvenecK. 2020. Climate change is likely to alter future wolf – moose – for- est interactions at Isle Royale National Park, United States. Frontiers in Ecology and Evolution 8:543915. de Pellegrin llorenTe, i., J. FAusKee, s. Burns and d. decKArd. 2020. Minnesota’s Forest Resources 2020. Minnesota Department of Natural Resources, Division of Forestry, St. Paul, Minnesota, USA. delgiudice, g.d., W.J. severud, T.r. oBermoller and B.d. smiTh. 2019. Winter nutritional restriction and decline of moose in northeastern Minnesota, win- ters 2013-2019. Saint Paul, USA: Minnesota Department of Natural Resources. Pages 196-212 in L. Cornicelli, M. Carstensen, M. Larsen, N. Davros, and B. Davis, editors. Summaries TWIG PRODUCTION PATTERNS ALCES VOL. 60, 2024 104 of wildlife research findings, 2018. Minnesota Department of Natural Resources, Wildlife Populations and Research Unit, St. Paul, Minnesota, USA. dreW, m.l., W.m. sAmuel, g.m. luKiWsKi and J.n. WillmAn. 1985. An evaluation of burning for control of winter ticks, Dermacentor albipictus, in central Alberta. Journal of Wildlife Diseases 21:313–315. ellingWood, d.d., P.J. PeKins, h. Jones and A.r. musAnTe. 2020. Evaluating moose Alces alces population response to infes- tation level of winter ticks Dermacentor albipictus. Wildlife Biology 2020. FArmer, r.e. 1962. Aspen root sucker for- mation and apical dominance. Forest Science 8:403–410. FelTon, A.m., e. holmsTröm, J. mAlmsTen, A. FelTon, J.P.g.m. cromsigT, l. edenius, g. ericsson, F. Widemo and h.K. WAm. 2020. Varied diets, including broadleaved forage, are important for a large herbivore species inhabiting highly modified landscapes. Scientific Reports 10:1904. Fisher, J.T. and l. WilKinson. 2005. The response of mammals to forest fire and timber harvest in the North American boreal forest. Mammal Review 35:51–81. gArWood, T.J., s.A. moore, n.m. FounTAin- Jones, P.A. lArsen and T.m. WolF. 2023. Species in the feces: DNA metabarcoding to detect potential gas- tropod hosts of Parelaphostrongylus tenuis consumed by moose (Alces alces). Journal of Wildlife Diseases 59:640–50. gArWood, T.J., W.J. severud, s.K. Windels, A.WünschmAnn, e.J. isAAc, A.g. Armien, s.A. moore and T.m. WolF. 2025. Predation vs. parasitism: A case study of Indigenous co-stew- ardship and science co-production to measure temporal shifts in moose mor- tality on ancestral lands of the Grand Portage Ojibwe. Ecology and Evolution:In Press grigAl, d.F. and l.F. ohmAnn. 1977. Biomass Estimation for some Shrubs from Northeastern Minnesota. U.S. Department of Agriculture Forest Service, Northern Research Station, Research Note NC-226. guPTA, r.d. and d. Kundu. 2001. Exponentiated exponential family: An alternative to gamma and Weibull distri- butions. Biometrical Journal 43:117–30. hArPer, g.J. 2008. Quantifying branch, crown and bole development in Populus tremuloides Michx. from north-eastern British Columbia. Forest Ecology and Management 255: 2286–2296. hAyes, F.P., J.J. millsPAugh, e.J. BergmAn, r.m. cAllAWAy and c.J. BishoP. 2022. Effects of willow nutrition and morphol- ogy on calving success of moose. The Journal of Wildlife Management 86:e22175. hoAglAnd, s.J. and s. AlBerT. 2023. Wildlife Stewardship on Tribal Lands: Our Place Is in Our Soul. Johns Hopkins University Press, Baltimore, MD. USA. hodgson, A.l. 2010. Temporal changes in spatial patterns of moose browse, causes and consequences. Ph.D. Thesis. University of Minnesota. 156 pp. homer, c., J. deWiTz, l. yAng, s. Jin, P. dAnielson, g. xiAn, J. coulsTon, n. herold, J. WicKhAm, and K. megoWn. 2015. Completion of the 2011 National Land Cover Database for the Conterminous United States – Representing a Decade of Land Cover Change Information. Photogrammetric Engineering and Remote Sensing 81:345–354. hoThorn, T., F. BreTz, P. WesTFAll, r.m. heiBerger, A. schueTzenmeisTer and s. scheiBe. 2023. multcomp: Simultaneous Inference in General Parametric Models. (version 1.4-25). https://cran.r-project. org/web/packages/multcomp/index. html. https://cran.r-project.org/web/packages/multcomp/index.html https://cran.r-project.org/web/packages/multcomp/index.html https://cran.r-project.org/web/packages/multcomp/index.html ALCES VOL. 60, 2024 TWIG PRODUCTION PATTERNS 105 hoy, s.r., l.m. vuceTich, r.o. PeTerson and J.A. vuceTich. 2021. Winter tick burdens for moose are positively associ- ated with warmer summers and higher predation rates. Frontiers in Ecology and Evolution 9:758374. hoy, s. r., J. s. ForBey, d. P. melody, l. m. vuceTich, r. o. PeTerson, K.B. KoiTzsch, l. o. KoiTzsch, A. l. von duyKe, J. J. henderson, g. l. PAriKh, and J.A. vuceTich. 2021. The nutritional condition of moose co-varies with cli- mate, but not with density, predation risk or diet composition. Oikos, 2022:e08498. irWin, L.L. 1985. Foods of moose, Alces alces, and white-tailed deer, Odocoileus virginianus, on a burn in boreal forest. Canadian Field Naturalist 99:240-245. Johnson, c.J. and r.v. reA. 2024. Response of moose to forest harvest and manage- ment: a literature review. Canadian Journal of Forest Research 54:366–88. Johnson, l.B., e. lArson, K.g. gill and J. sAvAge. 2022. Blending tree-ring fire- scar records and Indigenous memory in northern Minnesota, USA. Past Global Changes Magazine 30:36-37. KAssAmBArA, A. 2023. rstatix: Pipe-Friendly Framework for Basic Statistical Tests. (version 0.7.2). https://cran.r-project. org/web/packages/rstatix/index.html. Kimmerer, r.W. and F.K. lAKe. 2001. The role of Indigenous burning in land man- agement. Journal of Forestry 99:36–41. KrAFT, g.J., d.J. mechenich, c. mechenich, J. cooK and s.m. seiler. 2010. Assessment of natural resource conditions: Isle Royale National Park. Natural Resource Report NPS/NRPC/WRD/NRR. Department of the Interior, National Park Service, Fort Collins, CO. USA. KrAFT, g.J., d.J. mechenich, c. mechenich, m.d. WATerhouse, J. mcnelly, J. dimicK and J. cooK. 2014. Natural resource condition assessment: Grand Portage National Monument (revised July 2014). Natural Resource Report NPS/GRPO/NRR. Department of the Interior, National Park Service, Fort Collins, CO. USA. KreFTing, l.W., h.l. hAnsen and m.h. sTenlund. 1956. Stimulating regrowth of mountain maple for deer browse by herbicides, cutting, and fire. Journal of Wildlife Management 20: 434–41. lArson, e.r., K.F. KiPFmueller and l.B. Johnson. 2021. People, fire, and pine: linking human agency and landscape in the Boundary Waters Canoe Area Wilderness and beyond. Annals of the American Association of Geographers 111:1–25. lAuTenschlAger, R.A. 1992. Effects of conifer release with herbicides on moose: browse production, habitat use, and residues in meat. Alces 28: 215–22. leAK, W.B., d.s. solomon and s.m. FiliP. 1969. A silvicultural guide for northern hardwoods in the northeast. Upper Darby, PA: U.S. Department of Agriculture Forest Service, Northeastern Forest Experiment Station Research Paper NE-143. lenArz, M.S. 2007. 2007 aerial moose sur- vey. Saint Paul, USA: Minnesota Department of Natural Resources. lindroTh, R.L. 2001. Adaptations of quaking aspen for defense against damage by her- bivores and related environmental agents. Pages 273-284 in Shepperd, W. D., D. Binkley, D. L. Bartos, T. J. Stohlgren, J. Thomas, and L. G. Eskew, compilers. Proceedings: sustaining aspen in western landscapes: 2000, June13-15, Grand Junction, CO. Proceedings RMRS-P-18. Rocky Mountain Research Station, U.S. Department of Agriculture Forest Service, Fort Collins, CO. lorAnger, A.J., T.n. BAiley and W.W. lArned. 1991. Effects of forest succes- sion after fire in moose wintering habi- tats on the Kenai Peninsula, Alaska. Alces 27:100–109. lusK, c.h., P.B. reich, r.A. monTgomery, d.d. AcKerly and J. cAvender-BAres. 2008. Why are evergreen leaves so https://cran.r-project.org/web/packages/rstatix/index.html https://cran.r-project.org/web/packages/rstatix/index.html TWIG PRODUCTION PATTERNS ALCES VOL. 60, 2024 106 contrary about shade? Trends in Ecology and Evolution 23:299–303. mAngiAFico, S. 2024. rcompanion: functions to support extension education program evaluation. (version 2.4.35). https:// cran.r-project.org/web/packages/rcom- panion/index.html. mAnly, B.F., l. mcdonAld, d. ThomAs, T.l. mcdonAld and W.P. ericKson. 2002. Resource Selection by Animals: Statistical Design and Analysis for Field Studies. Second Edition. Kluwer Academic Publishers, New York, N.Y. USA. mAsTenBrooK, B. and h. cumming. 1989. Use of residual strips of timber by moose within cutovers in northwestern Ontario. Alces 25:146–55. mAThisen, K.M., J.M. Milner and C. Skarpe. 2017. Moose–tree interactions: rebrows- ing is common across tree species. BMC Ecology 17:12. mcnicol, J.g. and F.F. gilBerT. 1980. Late winter use of upland cutovers by moose. Journal of Wildlife Management 44:363–71. minnesoTA dnr ForesTry. 2016. Regeneration Monitoring - Procedures and Standards. St. Paul, MN. USA. 24 pp. miquelle, d.g. and P.A. JordAn. 1979. The importance of diversity in the diet of moose. Alces 15:54-79. moore, s.A., W.J. severud, T.m. WolF, K. PelicAn, J. BAuerKemPer, m. cArsTensen and s.K. Windels. 2024. Indigenous co-stewardship of North American moose: recommendations and a vision for a restoration framework. Journal of Wildlife Management 88:e22623. mummA, m.A., A.r. BevingTon, s. mArshAll and m.P. gillinghAm. 2024. Delineating wildfire burns and regrowth using satel- lite imagery to assess moose (Alces alces) spatial responses to burns. Ecosphere 15:e4793. nAdArAJAh, S. 2011. The exponentiated exponential distribution: a survey. Advances in Statistical Analysis 95:219–51. nAsh, J.C. 2014. On best practice optimiza- tion methods in R. Journal of Statistical Software 60: 1–14. nATure conservAncy. 1999. Classification of the Vegetation of Isle Royale National Park. USGS-NPS Vegetation Mapping Program. https://irma.nps.gov/DataStore/ DownloadFile/424670. neKolA, J.C. 2002. Effects of fire manage- ment on the richness and abundance of central North American grassland land snail faunas. Animal Biodiversity and Conservation 25:53-66. oliveirA-sAnTos, l.g.r., s.A. moore, W.J. severud, J.d. ForesTer, e.J. isAAc, y. chenAux-iBrAhim, T. gArWood, l.e. escoBAr and T.m. WolF. 2021. Spatial compartmentalization: A nonlethal pred- ator mechanism to reduce parasite trans- mission between prey species. Science Advances 7:eabj5944. PAinTer, l.e., r.l. BeschTA, e.J. lArsen and W.J. riPPle. 2018. Aspen recruit- ment in the Yellowstone region linked to reduced herbivory after large carnivore restoration. Ecosphere 9: e02376. PArdoe, I. 2020. Applied Regression Modeling. John Wiley and Sons, Hoboken, NJ. USA. PAriKh, g.l., J.s. ForBey, B. roBB, r.o. PeTerson, l.m. vuceTich and J.A. vuceTich. 2016. The influence of plant defensive chemicals, diet composition, and winter severity on the nutritional condition of a free-ranging, generalist herbivore. Oikos 126:196-203. PAsTor, J. 2016. What Should a Clever Moose Eat?: Natural History, Ecology, and the North Woods. Island Press, Washington, D.C. USA. PeeK, J.m., d.l. urich and r.J. mAcKie. 1976. Moose habitat selection and rela- tionships to forest management in north- eastern Minnesota. Wildlife Monographs 48:3-65. PerAlA, D.A. 1977. Manager’s handbook for aspen in the north-central states. U.S. Department of Agriculture Forest https://cran.r-project.org/web/packages/rcompanion/index.html https://cran.r-project.org/web/packages/rcompanion/index.html https://cran.r-project.org/web/packages/rcompanion/index.html https://irma.nps.gov/DataStore/DownloadFile/424670 https://irma.nps.gov/DataStore/DownloadFile/424670 ALCES VOL. 60, 2024 TWIG PRODUCTION PATTERNS 107 Service, North Central Forest Experiment Station General Technical Report NC-36. PeTerson, s., d. KrAmer, J. hursT and J. FrAir. 2020. Browse selection by moose in the Adirondack Park, New York. Alces 56:107–26. PeTerson, s., d. KrAmer, J. hursT, d. sPAlinger and J. FrAir. 2022. Forage and habitat limitations for moose in the Adirondack Park, New York. Alces 58:1-30. PorTingA, r.l.W. and r.A. moen. 2015. A novel method of performing moose browse surveys. Alces 51:107–122. PueTTmAnn, K.J. and P.B. reich. 1995. The differential sensitivity of red pine and quaking aspen to competition. Canadian Journal of Forest Research 25: 1731–1737. r core TeAm. 2024. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/. renecKer, l.A. and r.J. hudson. 1988. Seasonal quality of forages used by moose in the aspen-dominated boreal forest, central Alberta. Holarctic Ecology 11:111-118. renecKer, l.A. and c.c. schWArTz. 2007. Food Habits and Feeding Behavior. in C. C. Schwartz, A. W. Franzmann and R. E. McCabe, editors. Ecology and Management of the North American Moose. Second Edition. University Press of Colorado, Denver, CO. USA. ricArd, J.-g. and g.J. douceT. 1999. Winter use of powerline rights-of-way by moose (Alces alces). Alces 35:31–40. risenhoover, K.L. 1987. Winter foraging strategies of moose in subarctic and boreal forest habitats. Ph.D. Thesis, Michigan Technological University, MI, USA. roBBins, c.T., T.A. hAnley, A.e. hAgermAn, o. hJelJord, d.l. BAKer, c.c. schWArTz and W.W. mAuTz. 1987. Role of tannins in defending plants against ruminants: Reduction in protein availability. Ecology 68:98–107. rouTledge, r.g. and J. roese. 2004. Moose winter diet selection in central Ontario. Alces 40:95–101. sAFFord, L.O. 1983. Silvicultural guide for paper birch in the northeast (revised). Broomall, PA: U.S. Department of Agriculture Forest Service, Northeastern Forest Experiment Station Research Paper NE-535. sAnders, s. and J. KirschBAum. 2023. Woody species response to altered her- bivore pressure at Isle Royale National Park. Ecosphere 14:e4623. scheiner, s.m., T.l. shAriK, m.r. roBerTs and r.v. KoPPle. 1988. Tree density and modes of tree recruitment in a Michigan pine-hardwood forest after clear-cutting and burning. Canadian Field Naturalist 102:634–638. schrAge, M. 2023. 2023 Moose Habitat Survey. Fond du Lac Resource Management Division. Cloquet, MN. USA. schremPP, T.v., J.l. rAchloW, T.r. Johnson, l.A. shiPley, r.A. long, J.l. Aycrigg and m.A. hurley. 2019. Linking forest management to moose population trends: The role of the nutritional land- scape. PLoS ONE 14:e0219128. severud, W.J., g.d. delgiudice, T.r. oBermoller, T.A. enrighT, r.g. WrighT and J.d. ForesTer. 2015. Using GPS col- lars to determine parturition and cause-specific mortality of moose calves. Wildlife Society Bulletin 39:616–625. severud, W.J., T.r. oBermoller, g.d. delgiudice and J.r. FieBerg. 2019. Survival and cause-specific mortality of moose calves in northeastern Minnesota. Journal of Wildlife Management 83:1131–42. severud, W.J., s.s. Berg, c.A. ernsT, g.d. delgiudice, s.A. moore, s.K. Windels, r.A. moen, e.J. isAAc and T.m. WolF. 2022. Statistical population reconstruc- tion of moose (Alces alces) in northeast- ern Minnesota using integrated population models. PLoS ONE 17: e0270615. https://www.R-project.org/ TWIG PRODUCTION PATTERNS ALCES VOL. 60, 2024 108 severud, W.J., m.P. giguere, T. WAlTers, T.J. gArWood, K. TeAger, K.m. mArcheTTo, l.g.r. oliveirA-sAnTos, s.A. moore and T.m. WolF. 2023. Terrestrial gastropod species-specific responses to forest management: Implications for Parelaphostrongylus tenuis transmission to moose. Forest Ecology and Management 529:120717 smiTh, W.R. 2008. Trees and Shrubs of Minnesota. First Edition. University of Minnesota Press, Minneapolis, MN. USA. sTreeT, g.m., J. FieBerg, A.r. rodgers, m. cArsTensen, r. moen, s.A. moore, s.K. Windels and J.d. ForesTer. 2016. Habitat functional response mitigates reduced foraging opportunity: implica- tions for animal fitness and space use. Landscape Ecology 31:1939–53. sTulTs, m., s. PeTersen, J. Bell, W. BAule, e. nAsser, e. giBBons and m. FougerAT. 2016. Climate Change Vulnerability Assessment and Adaptation Plan: 1854 Ceded Territory Including the Bois Forte, Fond du Lac, and Grand Portage Reservations. 1854 Treaty Authority, Duluth, MN. USA. TerrA-Berns, m.h. 1993. Quantification and comparison of terpene concentra- tions in various balsam fir growth forms and foliage ages, and a simulation of moose browsing on balsam fir trees at Isle Royale. Unpublished M.S. thesis, Texas A&M University. 46 pp. TheriAulT, M.K. 2006. Moose to Moccasins: The Story of Ka Kita Wa Pa No Kwe. 2 edition. Dundurn Press, Toronto, Ontario. Canada. u.s. nATionAl PArK service. 2019. Weekly List 20190201 - National Register of Historic Places. https://www.nps.gov/ sub jec t s /na t iona l r eg i s t e r /week- ly-list-20190201.htm. vuKelich, J. 2023. The Seven Generations and The Seven Grandfather Teachings. Self-published. https://www.jamesvuke- lich.com/books. WAlTers, m.B. and P.B. reich. 1999. Low- light carbon balance and shade tolerance in the seedlings of woody plants: do winter deciduous and broad-leaved ever- green species differ? The New Phytologist 143:143–54. WAm, h.K., A. m. FelTon, c. sTolTer, l. nyBAKKen, and o. hJelJord. 2018. Moose selecting for specific nutritional composition of birch places limits on food acceptability. Ecology and Evolution 8: 1117–1130. Wehr, n.h., s.A. moore, e.J. isAAc, K.F. Kellner, J.J. millsPAugh and J.l. BelAnT. 2024. Moose and white-tailed deer mortal- ity peaks in fall and late winter. Journal of Wildlife Management 88: e22580. Wisconsin dnr. 2021. Chapter 21: Natural Regeneration, Appendix A: Forest Regeneration Metric, Pages 32-35 in Wisconsin Silviculture Handbook. Wisconsin DNR, Madison, WI, USA. WolF, T.m., y.m. chenAux-iBrAhim, e.J. isAAc and s.A. moore. 2021. Neonate health and calf mortality in a declining population of North American moose (Alces alces americanus). Journal of Wildlife Diseases 57:40-50 WünschmAnn, A., A.g. Armien, e. BuTler, m. schrAge, B. sTromBerg, J.B. Bender, A.m. FirshmAn and m. cArsTensen. 2015. Necropsy findings in 62 opportu- nistically collected free-ranging moose (Alces alces) from Minnesota, USA (2003-2013). Journal of Wildlife Diseases 51:157–165. zenner, e.K. and K.J. PueTTmAnn. 2008. Contrasting release approaches for a mixed paper birch (Betula papyrifera)– Quaking aspen (Populus tremuloides) stand. Northern Journal of Applied Forestry 25: 124–32. zuur, A., e.n. ieno, n. WAlKer, A.A. sAveliev and g.m. smiTh. 2009. Mixed effects models and extensions in ecology with R. Springer Science and Business Media, New York, NY. USA. https://www.nps.gov/subjects/nationalregister/weekly-list-20190201.htm https://www.nps.gov/subjects/nationalregister/weekly-list-20190201.htm https://www.nps.gov/subjects/nationalregister/weekly-list-20190201.htm https://www.jamesvukelich.com/books https://www.jamesvukelich.com/books