PATTERN AND PROCESS OF TREE REGENERATION AND RECRUITMENT IN MANAGED NORTHERN HARDWOOD FORESTS By Catherine Rose Henry A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Forestry—Doctor of Philosophy Ecology, Evolutionary Biology and Behavior—Dual Major 2022 ABSTRACT PATTERN AND PROCESS OF TREE REGENERATION AND RECRUITMENT IN MANAGED NORTHERN HARDWOOD FORESTS By Catherine Rose Henry For managed forests which rely on natural tree regeneration for canopy recruitment, abundance and composition of tree regeneration portend future forest structure and diversity. For northern hardwood forests, a geographically widespread forest type in North America, typical single-tree selection (STS) management relies on natural regeneration to promote new cohorts of canopy trees. Harvesting dispersed, select trees every 15 – 20 years, STS generates low light levels intended to promote sugar maple (Acer saccharum) and other shade-tolerant tree species in an uneven-aged system. However, following 60 + years of STS implementation in the Great Lakes region, concerning regeneration trends have emerged, namely low densities of sugar maple and low tree species diversity. Additionally, few studies have analyzed age structure under this system to assess its past efficacy in generating uneven-aged forests. The research presented here characterizes regeneration and recruitment outcomes of STS, analyzing data from a uniquely detailed and geographically widespread research project of 141 northern hardwood stands across northern Michigan. Given the silvicultural focus on regenerating sugar maple, the first two chapters focus on management outcomes for this key species. First, a flexible Bayesian hierarchical model offers insight on patterns of sugar maple regeneration for key size classes as a function of plot and stand level predictors. Our results indicate that sugar maple regeneration is sparse to absent, particularly for size classes actively browsed by deer and recently escaped from the deer browsing zone. The second analysis characterizes age structure for a subset of 51 stands, drawing on 1499 sugar maple trees > 5 cm diameter sampled via basal discs from recently harvested stumps; this analysis provides insight to past patterns of recruitment and establishment. The results suggest little evidence of sugar maple seedling regeneration and canopy ingrowth over the past 60 + years of STS management; instead, stands have highly suppressed saplings plus aging poletimber and sawtimber classes, which are at or quickly approaching economic maturity. Given declines in sugar maple dominance as evidenced by the first two research analyses, the third analysis assesses stand-level tree species diversity and individual species abundance as a function of landscape predictors and size class to shed light on projected future canopy composition. On average, there are approximately three effective common species for seedlings, saplings, and canopy stems at the stand level, and species less desirable for management are occupying growing space in the sapling layer. Together, these results indicate that STS has been unsuccessful in regenerating or recruiting sugar maple over the past 60 + years, and stands are characterized by a paucity of tree species. Our results support several potential alternative management strategies, including decreasing basal area via more intense harvests, prohibiting deer browsing via natural browsing barriers, or introducing greater diversity of tree species via direct seeding or planting. These results should be considered to improve current management of northern hardwood forests in the Great Lakes region. To all who supported me along the way, especially my partner and my parents iv ACKNOWLEDGEMENTS First and foremost, I would like to thank my advisor, Dr. Mike Walters, for bringing me on to the Big Hardwoods project. Your continued support, guidance, and enthusiasm made all the difference. Thank you to my committee members, Dr. Andrew Finley, Dr. Richard Kobe, and Dr. Gary Roloff, for invaluable feedback and guidance on all things statistical, biological, and ecological. To my lab mates, Lucas Elenitsky and Evan Farinosi, I am honored to have collaborated and learned with you. Thank you to our wonderful field technicians, particularly to those who helped me cut and haul 1500+ tree discs out of the woods. My work would not have been possible without funding and support from the Michigan Department of Natural Resources. To our collaborators, especially Mike Donovan, Jason Hartman, and Tim Greco, thank you for making this project possible. My graduate studies were also funded by the Michigan State University Plant Science Fellowship, the Department of Forestry Carol C. Gustafson and Gary S. Murphy Endowed Scholarship, the Department of Forestry Charles Guenther Memorial scholarship, and the Michigan State University Graduate School Writing Fellowship, for which I am most grateful. Thank you to my Michigan friends who became like family, you kept me grounded. To my parents and sisters, thank you for the love and support from the very beginning. And to my partner, without whom I would not have made it this far, thank you for everything. v TABLE OF CONTENTS LIST OF TABLES ..................................................................................................................... viii LIST OF FIGURES ................................................................................................................... xiii CHAPTER 1 INTRODUCTION ................................................................................................. 1 CHAPTER 2 COMPLEX DRIVERS OF SUGAR MAPLE (ACER SACCHARUM) REGENERATION REVEAL CHALLENGES TO LONG-TERM SUSTAINABILITY OF MANAGED NORTHERN HARDWOOD FORESTS .......................................................................................................................... 7 2.1 Abstract.......................................................................................................................... 7 2.2 Introduction ................................................................................................................... 8 2.2.1 Management and regeneration dynamics ............................................................... 8 2.2.2 Local and landscape factors affecting regeneration ............................................ 11 2.3 Methods ........................................................................................................................ 14 2.3.1 Study area ............................................................................................................. 14 2.3.2 Field methods ........................................................................................................ 16 2.3.3 Sugar maple regeneration density modeling framework ...................................... 17 2.4 Results .......................................................................................................................... 22 2.4.1 Sugar maple regeneration characteristics ............................................................ 22 2.4.2 Stand characteristic covariation ........................................................................... 26 2.4.3 Sugar maple regeneration models ........................................................................ 26 2.5 Discussion..................................................................................................................... 33 2.6 Implications for Management .................................................................................... 36 2.7 Acknowledgements ..................................................................................................... 40 CHAPTER 3 SUGAR MAPLE AGE STRUCTURE IN NORTHERN HARDWOOD STANDS MANAGED BY SELECTION SILVICULTURE ...................................... 41 3.1 Abstract........................................................................................................................ 41 3.2 Introduction ................................................................................................................. 42 3.3 Methods ........................................................................................................................ 50 3.3.1 Sites and study area .............................................................................................. 50 3.3.2 Field methods ........................................................................................................ 51 3.3.3 Lab methods .......................................................................................................... 52 3.3.4 Data analysis ........................................................................................................ 53 3.3.4.1 Estimating age and diameter at breast height .................................................. 53 3.3.4.2 Age associations with tree size characteristics (Hypothesis 1) ........................ 54 3.3.4.3 Age cohort analysis (Hypothesis 2) .................................................................. 55 3.3.4.4 Comparison with theoretical age structure (Hypothesis 3) .............................. 57 3.3.4.5 Stand and landscape factors associated stand age structure (Hypothesis 4) ... 58 3.4 Results .......................................................................................................................... 60 3.4.1 Data overview ....................................................................................................... 60 3.4.2 Associations between age and size characteristics (Hypothesis 1) ...................... 61 vi 3.4.3 Age cohort results (Hypothesis 2) ......................................................................... 64 3.4.4 Age structure deviation from expectations (Hypothesis 3) ................................... 66 3.4.5 Factors associated with stand age structure (Hypothesis 4) ................................ 67 3.5 Discussion..................................................................................................................... 70 3.5.1 Overview ............................................................................................................... 70 3.5.2 Comparison to other studies ................................................................................. 73 3.5.3 Management implications ..................................................................................... 75 3.5.4 Caveats and future directions ............................................................................... 76 CHAPTER 4 PATTERNS AMONG SPECIES FROM SEEDLING TO CANOPY PORTEND FUTURE COMPOSITIONAL SHIFTS AND LOW RESILIENCE IN MANAGED NORTHERN HARDWOOD FORESTS ................. 78 4.1 Abstract........................................................................................................................ 78 4.2 Introduction ................................................................................................................. 79 4.3 Methods ........................................................................................................................ 85 4.3.1 Study area and data collection ............................................................................. 85 4.3.2 Species diversity analyses ..................................................................................... 89 4.3.3 Generalized linear modeling for species abundance ............................................ 92 4.3.4 Relative abundance as a function of size class ..................................................... 92 4.4 Results .......................................................................................................................... 94 4.4.1 Overview of stand characteristics ......................................................................... 94 4.4.2 Species diversity measures predicted by stand and landscape factors ................. 96 4.4.3 Individual tree species density predicted by stand and landscape factors ........... 99 4.4.4 Species relative abundance across size classes .................................................. 101 4.5 Discussion................................................................................................................... 105 4.5.1 Overview ............................................................................................................. 105 4.5.2 Management implications ................................................................................... 111 4.5.3 Caveats ................................................................................................................ 112 CHAPTER 5 CONCLUSIONS................................................................................................ 114 APPENDICES ........................................................................................................................... 118 APPENDIX A Details on the deer model ............................................................................. 119 APPENDIX B Additional information on site quality categorizations ................................ 125 APPENDIX C Sugar maple model: parameter and model selection .................................... 130 APPENDIX D Details on species diversity analysis ............................................................ 138 APPENDIX E Model selection for stand-level tree species abundance ............................... 148 WORKS CITED........................................................................................................................ 150 vii LIST OF TABLES Table 2.1. For our multilevel models of the plot-level count of sugar maple regeneration (seedlings, small saplings, and large saplings), we included predictor variables at both plot and stand levels. The table lists a brief description, the abbreviation used, the mean, and the range of parameter values........................................................... 19 Table 2.2. The mean, median, and range of sugar maple regeneration densities for seedlings, small saplings, and large saplings at the plot and stand level. Seedlings are 0-50 cm tall, small saplings are 50 cm - 137 cm, and large saplings are 137 cm tall or greater and up to 5 cm DBH. Densities are presented in stems per plot (2 m2 plot for seedlings, 12.5 m2 plot for small and large saplings). Plot occurrences indicates the number of surveyed plots which contained at least one sugar maple. Stand occurrences indicates the number of sites which contained at least one sugar maple stem in any of the plots surveyed in that stand. The percent relative abundance refers to the percentage of sugar maple regeneration stems relative to the total density of stems in that size class; for plot percentage relative abundance, plots with zero total stems in a size class were omitted. ................................................... 24 Table 3.1. Summary of published sugar maple age-diameter analyses, including source, location, stand history, sampled DBH range, sampling height, the number of stands (nstands), the number of trees sampled (nsamples), the form of the age-diameter relationship, estimated age of 10 cm DBH stems, and the reported R2. Stand history codes include old growth (OG), selection management (SEL), partial cuttings (PART), variable histories, including SEL and PART (VAR), and unmanaged (UM). For sampling height, S refers to stump height (unspecified). Age-diameter relationships are linear (L), quadratic (Q), or free-hand curve (C, which is not a statistical model and therefore does not have an R2 (N/A)). Tubbs (1977a) contained two snapshots of the same stand, before and after harvest treatment, and so is included twice. Harmala (2021) reported stand models separately; for brevity, we have summarized the unmanaged control from the harvest-managed stands, but additional details can be found in Harmala (2021). Information not reported is marked with a *, and information not applicable is marked N/A. ...................................................................................................................... 44 Table 3.2. For predictors of AGE10, AGE20, and AGE30, mean and range of for continuous variables and number for categorical site quality. Site quality 1 is lowest and site quality 4 is highest. .................................................................................. 59 Table 3.3. Plot-level diameter composition of all living stems > 10 cm DBH for 51 plots, each representative of stand, in which age discs were sampled. ...................................... 61 Table 3.4. Model parameter estimates (E 50%) and 95% credible intervals (E 2.5%, E 97.5%) for factors predicting age of sugar maple stems at 10, 20, and 30 cm DBH. Significant predictors (p < 0.05) are bolded and italicized. Site quality 1 (lowest viii site quality) is included in the intercept, and site quality 4 is the highest site quality. T2 represents the variance (equal to the squared standard deviation). ................. 69 Table 4.1. Summary statistics for tree species regeneration at the stand level when present in a given size class, including mean, maximum (stems ha-1), and number of stand occurrences (n, out of 141) by size class, plus overall number of site occurrences n (across all size classes). Summary statistics were calculated for the subset of stands which that species was present on, in a given size class. Seedlings are 0 – 137 cm tall, small saplings are > 137 cm tall and < 10 cm DBH, and canopy trees are > 10 cm DBH. Species labeled with ** were modeled for individual abundance and relative abundance by size class; species labeled with * were modeled only for relative abundance by size class. Species without any * were included only for general species diversity calculations and stand summaries. ............... 87 Table 4.2. Summary statistics for predictor variables included seedling and sapling models for diversity and individual species distribution. Variables include estimated deer use (DEER), average basal area (m2ha-1, BA), standard deviation of basal area (BASD), seed source as estimated by the density of conspecific stems > 10 cm DBH (SeedS), estimated asymptotic species richness of stems > 10 cm DBH, estimated effective number of common species (Hill q=1) of stems > 10 cm DBH, % coverage of shrubs (shrub), % coverage of bare mineral soil (BMS), hardwood litter (HWL), and coarse woody debris (CWD). † SeedS was only included in the individual species models. * SHRUB, BMS, HWL, and CWD were only used as predictors for the seedling models....................................................... 91 Table 4.3. Parameter estimations for models of species richness (SR) and diversity (CS) for seedlings (0 – 137 cm tall), saplings (> 137 cm tall and < 10 cm DBH), and canopy/mature trees (> 10 cm DBH), as predicted by total canopy/mature basal area, standard deviation of canopy/mature basal area, and deer use, plus shrub coverage and percent hardwood litter coverage (HWL), coarse woody debris (CWD), and bare mineral soil (BMS) for the seedling models. Model fit parameters (Intercept, Beta and Shape) are also provided: Beta characterizes dispersion in a negative binomial distribution (SR model only) and shape influences the shape of the gamma distribution (CS model only). Significant associations (95% CI) are bolded. .................................................................................... 97 Table 4.4. Summaries of stand-specific relative abundance rankings by size class for 18 species surveyed in > 10 % of our 141 stands. Scientific name and authorship can be looked up by species code in Table 4.1. Species are grouped by the shape of their average relative abundance size class distribution, and n is the number of sites that a species was present on, in any size class. Seedling, sapling, and canopy % indicates the percentage of stands in which that size class had the highest relative abundance (seedling 0 – 137 cm tall, saplings > 137 cm and < 10 cm DBH, and canopy/mature, abbreviated as canopy, > 10 cm DBH). 1st and 2nd ranking refer to two most common rankings of relative abundances by size class (n is number of stands). The nomenclature presents the size classes in decreasing rank order (i.e., the first size category listed had the highest relative abundance), ix with size classes which had zero stems omitted. For example, 21 out of 77 stands in which balsam fir (BF) was present had the relative abundance ranking “Can- Sap-Seed” (canopy is highest relative abundance); as a left-skewed species, it has greatest relative abundance in the canopy, on average. .................................................. 104 Table A.1. Results from the stage one deer use model are presented with their median value and a 95% Bayesian confidence interval. Sigma squared is the variance parameter for the spatial term, and effective spatial range represents the distance at which spatial autocorrelation becomes insignificant in the model, in kilometers. None of the climatic or landscape predictor variables were significantly different from zero, likely because the spatial term explained most of the variation. ................... 124 Table B.1. Bayesian multilevel model predicting sugar maple regeneration in Michigan included four ordinal categories of site quality (sensu (Burger & Kotar, 2003). Kotar site quality grouping, model abbreviation, Kotar types (see Table B.2) included in the grouping, and count of stands sampled by region. The full set of site qualities was only present in one of the regions. NLP = Northern Lower Peninsula, EUP = Eastern Upper Peninsula, WUP = Western Upper Peninsula. *Medium site quality was included in the intercept. ...................................................... 126 Table B.2. Descriptions of site quality types (Burger and Kotar 2003) used to model sugar maple regeneration. NLP = northern Lower Peninsula, WUP = western Upper Peninsula, and EUP = eastern Upper Peninsula................................................... 127 Table C.1. Adjusted R2 values for a linear models of stand-level log average sugar maple seedling counts predicted by different measures of canopy sugar maple density. Predictors included stand-level average stem count and basal area of sugar maple trees within 6 m2 radius plots, at varying size thresholds. The basal area of sugar maple trees > 25 cm had the highest adjusted R2 model and was used as the predictor of seed availability in model used to predict sugar maple regeneration. ......... 131 Table C.2. Results of the pairwise Wilcoxon tests using a Bonferroni correction (n=24). Values indicate probability of the distributions of the pair being significantly different. Values significant at adjusted p-values of 0.05 are highlighted. ..................... 132 Table C.3. Candidate models for sugar maple seedlings (SM1), small saplings (SM2), and large saplings (SM3) predicted by plot level densities of total small sapling stems (tot2), total large sapling stems (tot3), stems 5-10 cm DBH (sap), and shrub coverage (wc), plus the stand-level measurements of basal area of sugar maple stems > 25 cm DBH (SM10BA), total stand basal area > 10 cm (standBA), deer use (deer), and site quality (Nut1, Nut3, Nut4) which generate a stand-level intercept, a. The DIC values are based on 10,000 post burn-in samples from three MCMC chains using the RJAGS package (Plummer, 2018). For each size class, the model with the lowest DIC value is bolded. ............................................................. 136 Table D.1. Summary of stand (n = 51) classification under an herbaceous-indicator site- quality system (Burger & Kotar, 2003). Site quality 1 is lowest, 4 is highest. x Letters represent key indicator species: A = Acer saccharum, Aa = Aralia nudicaulis, Ar = Acer rubrum, As = Arisaema atrorubens, Ast = Aster marophyllus, Ca = Caulophyllum thalictroides, D = Dryopteris spinulosa, F = Fagus grandifolia, Hp = Hepatica [variant], Ly = Lycopodium annotinum, M = Maianthemum canadense, O = Osmorhiza claytoni [variant], P = Pinus strobus, Po = Polygonatum pubescens, Sm = Smilacina racemose [variant], T = Tsuga canadensis, V = Vaccinium angustifolium, Vb = Viburnum acerifolium, and [w] = Wisconsin variant............................................................................................................ 139 Table D.2. Summary of tree species (scientific and common name plus number of plot occurrences (out of 51)) from plot-level data, surveyed in 51 managed northern hardwood forests. ............................................................................................................ 140 Table D.3. AIC comparison between linear and Chapman-Richards (CR) models, by stand, plus the pseudo-R2 value of the selected model (lower AIC, selected model AIC is bolded). Stands are identified by region (EUP = eastern Upper Peninsula, WUP = western Upper Peninsula, and NLP = northern Lower Peninsula). ................... 141 Table D.4. Model parameter estimates, including mean (E), standard error (SE), and p- value (p-val), for AIC-selected Chapman Richards models. Model estimates significantly different from zero (p < 0.05) are bolded. The model equation for age of stem i as a function of DBH is: Agei ~ Asym*(1-e-b*DBHi )c , where the variable Asym is interpretable as the age at which model estimates plateau. Stands are identified by region (EUP = eastern Upper Peninsula, WUP = western Upper Peninsula, and NLP = northern Lower Peninsula).......................................................... 142 Table D.5. Model parameter estimates, including mean (E), standard error (SE), and p- value (p-val) for AIC-selected linear models. Model estimates significantly different from zero (p < 0.05) are bolded. The model (which is unique for each stand) for age of stem i a function of DBH is: Agei ~ Intercept + Slope*(DBHi ). Stands are uniquely numbered and identified by region (EUP = eastern Upper Peninsula, WUP = western Upper Peninsula, and NLP = northern Lower Peninsula). ....................................................................................................................... 144 Table D.6. Model summaries of fixed effects of crown class as a factor predicting sugar maple age (overtopped included in the intercept, intermediate CPIntermediate, codominant CPCodominant, and dominant CPDominant stems). Significant predictor estimates (p < 0.05) are bolded. ...................................................................................... 145 Table D.7. Model summary of the fixed effect variables crown percentage, DBH (cm), and their interaction as predictors of sugar maple age. Significant predictor estimates (p < 0.05) are bolded. ...................................................................................... 145 Table D.8. Model summary of a fixed effects of sugar maple stem age as a function of height. Significant predictor estimates (p < 0.05) are bolded. ........................................ 145 xi Table E.1. DIC comparison for negative binomial (NB) and zero-inflated negative binomial (ZINB) models. Lower DIC value is bolded. Species codes can be found in Table 4.1. Size class 1 refers to seedlings (0 – 137 cm tall) and size class 2 refers to saplings (> 137 cm tall and < 10 cm DBH). ..................................................... 149 xii LIST OF FIGURES Figure 2.1. Locations of 141 single-tree selection-managed northern hardwood stands sampled for sugar maple regeneration and recruitment across northern Michigan by site class. SQ1 stands are poor to poor/medium, SQ2 medium, SQ3 medium/rich, and SQ4 rich to very rich (Burger and Kotar, 2003). ................................. 16 Figure 2.2. Model estimated stand-level deer use, winter of 2017. Zero represents a stand with no estimated winter deer use (0 segments had fecal pellets), while 100 indicates a high use stand where all 103 transect segments (each segment ~6m x 6m) had fecal pellets. ........................................................................................................ 20 Figure 2.3. Average stand-level density of sugar maple seedlings (0 – 50 cm tall; A), small saplings (50 – 137 cm tall; B), and large saplings (>137 cm tall and <5 cm diameter breast height) based on data collected for 141 northern hardwood stands, summer 2017. .................................................................................................................... 25 Figure 2.4. Kendall’s tau correlation matrix of predictor variables used to model stand level sugar maple regeneration for Michigan (non-significant (p>0.05) correlations are labeled with an x, positive associations are blue, and negative associations are red with parallel lines). Variables include total stem density of small saplings (SAP1), total density of large saplings (SAP2), total density of trees 5 – 10 cm DBH (SUB), deer usage (DEER), site quality (SQ, treated as continuous from 1-4), shrub coverage (SHRUB), basal area of sugar maple trees > 25 cm (SM.BA), and total stand basal area (BA). ............................................................ 26 Figure 2.5. Model parameter estimates for predictors of sugar maple regeneration, with 95% Bayesian credible intervals. Variables include total stem density of small sugar maple saplings (SAP1), total density of large sugar maple saplings (SAP2), total density of trees 5 – 10 cm DBH (SUB), shrub coverage (SHRUB), deer usage (DEER), total stand basal area (BA), basal area of sugar maple trees > 25 cm (SM.BA), site quality (SQ1 is poor, SQ3 medium/rich, SQ4 rich to very rich, and SQ2 medium included in the intercept), and model intercept (MU), which represents average plot-level stem count, on the log scale, when all predictors are at average values. Quadratic terms (indicated by superscripted 2) included in the final models portrayed with dashed lines.......................................................................... 28 Figure 2.6. Model estimates for the relationship of sugar maple seedlings per plot (2 m2) with the plot-level densities of small saplings (SAP1), large saplings, (SAP2) and all sub-canopy trees 5 – 10 cm DBH (SUB); the stand-level basal area of sugar maple canopy trees larger than 25 cm DBH (SM.BA); and the plot-level density of shrubs (SHRUB). Mean (black line) and 95% Bayesian credible intervals (gray lines) are displayed, based on the average dispersion parameter estimate. ...................... 29 Figure 2.7. Model estimates for the relationship of sugar maple small saplings per plot (12.6 m2) with large saplings per plot (SAP2), subcanopy trees per plot (SUB), xiii deer use (DEER), and basal area of sugar maple canopy trees larger than 25cm DBH (SM.BA). Mean (black line) and 95% Bayesian credible intervals (gray lines) are displayed for the average dispersion parameter estimate. ................................. 30 Figure 2.8. Model estimates for the relationship of sugar maple large saplings per plot (12.6 m2) with deer use (DEER), basal area of sugar maple canopy trees larger than 25cm DBH (SM.BA), and shrub coverage (SHRUB). Mean (black line) and 95% Bayesian credible intervals (gray lines) are displayed for the average dispersion parameter estimate. .......................................................................................... 31 Figure 2.9. For each sugar maple regeneration class, boxplots of stand-level average density by site quality category (sensu Kotar and Burger 2003), northern Michigan, 2017. Site quality categories include poor (SQ1), medium (SQ2), medium to rich (SQ3), and rich to very rich (SQ4) in the model. Middle lines represent median values, the lower and upper hinge represent the first and third quartiles, and the whiskers represent values no further than 1.5 of the interquartile range (distance from first to third quartile); individual points are outliers. ...................... 32 Figure 2.10. Percentage of surveyed stands which would be considered stocked in small (50-137 cm tall) and large (137 cm tall to 5 cm DBH) sugar maple saplings, based on average stand densities, is plotted as a function of stocking criteria thresholds (100% stocking at a threshold of 0 stems ha-1, not shown). Red dashed lines provide an example of figure interpretation: at stocking criteria of 2500 stems ha-1 for small sugar maple saplings, 29% of surveyed stands are stocked, and at 1000 stems ha-1 for large sugar maple saplings, 35% of stands are stocked. ............................ 38 Figure 3.1. Stand-level diameter distribution of sugar maple, conifer, undesirable, and other desirable species in 5 cm diameter size classes (size classes are labeled with median diameter value) from 51 northern hardwood forests stands sampled in northern Michigan, USA. Undesirable includes American beech and ironwood, while other desirable refers to all species other than sugar maple and undesirable species (sensu Walters et al., 2022 but excluding ash as undesirable species). Black asterisks represent theoretical diameter class stem densities at full stocking, for all species > 5 cm DBH, according to Arbogast (1957).............................................. 60 Figure 3.2. Age versus diameter (cm) for 1499 sugar maple stems ≥ 5 cm in managed northern hardwood forests (black points). The Chapman Richards model fit to our data is shown in orange. The other lines represent relationships between age and diameter for sugar maple stems from other studies (for the range of diameter values in each study). Details on studies can be found in Table 3.1................................. 62 Figure 3.3. Model-predicted values of sugar maple age by diameter (for the stand-specific range of diameters) for 51 managed northern hardwood forests. Black lines are linear models and red lines are Chapman-Richards models. ............................................ 63 Figure 3.4. Histogram (5-year bins) and boxplot of median age for sampled discs in Arbogast’s sapling (5 – 11.4 cm DBH), poletimber (11.4 – 24.1 cm DBH), and xiv sawtimber classes (> 24.1 cm DBH), by stand, for 51 STS managed northern hardwood forests. Dashed lines emphasize 50, 100, and 150 years on the x axis. Boxplots show first, second (median), and third quartiles, plus whiskers which extend up to approximately three standard deviations from the mean (1.5 times the inter-quartile range); points represent outliers. Note, we did not sample the full width of Arbogast’s sapling class (stems 3.8 – 11.4 cm DBH). ....................................... 64 Figure 3.5. Histogram-identified sugar maple age cohorts by site (x-axis) for 51 managed northern hardwood stands (3.5A). Point size represents total stems per hectare measured within that cohort, and black points represent the age cohort capturing the dominant canopy. Stands are organized by number of identified cohorts within each stand and sorted (low to high) by median age of the youngest age cohort. Number of identified age cohorts (3.5B) and median age of the canopy cohort (3.5C) are represented spatially, plus a histogram of canopy ages (3.5C)........................ 66 Figure 3.6. Bar plot of average age density distribution (20-year age bins, labeled with median value) of sugar maple stems > 5 cm DBH, pooled across 51 managed NHF stands; black points represent age density of a fully stocked northern hardwood stand, following harvest, if entirely stocked with sugar maple (Figure 3.6A). Anticipated age distribution was determined by applying the age-diameter relationship from Tubbs (1977a) to the size class densities outlined by Arbogast (1957). Bar plot representing the percentage of stands individually considered stocked (e.g., exceed Arbogast’s extrapolated stocking value) for each age class expected under the anticipated age distribution (Figure 3.6B). ........................................ 67 Figure 3.7. Kernel density estimates of model-estimated values of stems 10 cm (purple solid line), 20 cm (green dashed line), and 30 cm (yellow dotted line) DBH are presented (n = 48 stands for AGE10, n = 50 for AGE20, and n = 51 for AGE30); these are interpretable as a smoothed version of a histogram, where the area under the curve of each line sums to one. We did not extract model estimates beyond the range of data, therefore estimates were not available for three stands lacking saplings, one of which additionally had no stems below 20 cm DBH. ............................ 68 Figure 4.1. Landscape-level (i.e., across 141 stands in Michigan) average relative abundance, by size class, for species which achieved at least 1% relative abundance in any size class; species include sugar maple (SM), American beech (AB), ironwood (IW), red maple (RM), balsam fir (BF), white ash (WA), cherry species (CH), basswood (BW), and yellow birch (YB). Scientific nomenclature and all species included in “Other” can be found in Table 4.1. Seedlings are 0 – 137 cm tall, saplings are > 137 cm tall and < 10 cm DBH, and canopy/mature are > 10 cm DBH. ................................................................................................................... 95 Figure 4.2. Stand-level estimated asymptotic species richness (SR) and effective number of common species (CS) for seedlings (0 –137 cm tall), saplings (> 137 cm tall and < 10 cm DBH), and canopy/mature trees (abbreviated canopy for brevity, > 10 cm DBH). Asymptotic estimates account for sampling effort and are the xv estimated true number of species present in the stand. Within each panel, box plots labeled with the same letter are not significantly different. ..................................... 96 Figure 4.3. Histograms for asymptotic species richness (SR) and number of common species, CS, by site quality and region. Site quality 1 = lowest quality and 4 = highest (Burger and Kotar 2003). Letters denote significant differences (pairwise comparison of estimated marginal means) among site qualities within individual panels; bars with the same letters are not significantly different. ..................................... 98 Figure 4.4. Parameter estimates (median with 0.05 – 0.95 confidence interval) for seedling (0 – 137 cm tall) and sapling (> 137 cm tall and < 10 cm DBH) densities with predictor variables by tree species: average stand canopy/mature basal area (BA), standard deviation of canopy/mature basal area (BASD), deer use (Deer), seed source as estimated by conspecific canopy/mature tree density (SeedS), shrub coverage (Shrub), percent bare mineral soil (BMS), percent hardwood litter (HWL), and percent coarse woody debris (CWD). CWD, HWL, BMS, and Shrub were only included in the seedling model. Median values of significant (p > 0.10) parameters are diamonds, outlined in black. Seedlings and saplings of AB, CH, IW, RM, and SM, plus BW seedlings were modeled with a negative binomial model; all others were modeled with a zero-inflated model. .......................................... 100 Figure 4.5. Relative abundance by species relative to stems of all species of the same size class for seedlings (0 –137 cm tall), saplings (> 137 tall and < 10 cm DBH), and canopy/mature (abbreviated canopy in figure, > 10 cm DBH) strata. Within a species plot, different letters indicate different significances in relative abundance between size classes (p-value < 0.05). Size class was not a significant predictor for American elm and white pine, so we did not calculate pairwise comparisons. Histogram color corresponds to their identified shape across size classes (legend on right). .......................................................................................................................... 102 Figure A.1. Design of the winter pellet survey transects. Points represent 25 sampling locations for vegetation surveys within the 12.14 ha square unit. Each transect (blue line) was 6 m wide, and cumulative transect length was 628 m. For analysis, transects were divided into spatially referenced segments ~ 6 m long. Occurrence of deer pellets were documented for each segment. ....................................................... 121 Figure C.1. Decision tree generated by ANOVA model recursive partitioning for total stand-level counts of sugar maple seedlings (C.1A), small saplings (C.1B), and large saplings (C.1C). Predictive inputs were unique identifiers for site quality (1st, 2nd, 3rd, 4th, with 1st being the highest quality and 4th being the lowest quality) by region (eastern Upper Peninsula - EUP, western Upper Peninsula - WUP, and northern Lower Peninsula - NLP); for example, 1st EUP indicates highest quality sites in the eastern Upper Peninsula. In the bottom boxes, count values indicate the predicted value (total count of sugar maple stems surveyed per stand) while percentage values indicate the percentage of observations (out of 141) in the node. ............................................................................................................. 133 xvi Figure D.1. Individual tree disc ages (n = 1499) by competitive class (O = overtopped, I = intermediate, C = codominant, and D=dominant). ......................................................... 146 Figure D.2. Map of model-averaged sugar maple stem age at 10 cm DBH for 51 managed northern hardwood forests. ............................................................................................. 146 Figure D.3. Maps of average sample age for sugar maple stems in 51 managed northern hardwood forests in three size classes: saplings (5 – 11.4 cm DBH), poletimber (11.4 – 23.9 cm DBH), and sawtimber (> 23.9 cm DBH). ............................................. 147 xvii CHAPTER 1 INTRODUCTION Forests comprise approximately 30% of the world’s total land area (Global Forest Resources Assessment, 2015) and provide a variety of cultural, economic, and ecological services (Brockerhoff et al., 2017). Globally, continued provisioning of these services relies on how forests respond to increasing human-mediated stressors, including increasing forest fire activity (Flannigan et al., 2000), invasive plants and pathogens (Dukes et al., 2009; Ramsfield et al., 2016), and frequency and severity of drought (Allen et al., 2010), plus interactions among factors (Dale et al., 2000) and with climate change (Dale et al., 2001). Given trees are sessile, relatively long-lived individuals with limited seed dispersal distances, migration rates have not matched climatic changes, nor are they projected to in the future (Iverson et al., 2004; Zhu et al., 2012), generating an adaptational lag (Aitken et al., 2008). Therefore, high forest biodiversity, which generally increases resiliency (Thompson et al., 2009), is key to future forest resiliency. Current patterns of low density and diversity of tree regeneration (Ramirez et al., 2018; Miller & McGill, 2019; Vickers et al., 2019) for North American forests are thus significant cause for concern. Covering over 50 million hectares of northeastern United States (Oswalt et al., 2014), northern hardwood forests are an important forest ecosystem. They are currently ecologically dominated by sugar maple (Schulte et al., 2007), which is also one of the most economically valuable species (Linehan & Jacobson, 2005; Duval et al., 2014). Northern hardwood forests can support a variety of species, including deciduous trees such as maples (Acer spp.), American beech (Fagus americana), oaks (Quercus spp.), basswood (Tilia americana), aspen (Populus spp.), ash (Fraxinus spp.), and elm (Ulmus spp.), in addition to conifers such as hemlock (Tsuga 1 canadensis), spruce (Picea spp.), and fir (Abies balsamifera) (Schulte et al., 2007). However, species richness and structural complexity have declined since European colonization (Schulte et al., 2007). Management of northern hardwood forests in the Great Lakes region has been dominated by a silvicultural system known as single-tree selection. Development of this system began in the 1920’s (Kern et al., 2014), in response to a preponderance of young, even-aged forests established following widespread clearcutting and often subsequent slash-fueled fires in the late 19th and early 20th centuries (Whitney, 1987; Dickmann & Leefers, 2016). Research and development of this silvicultural system (Eyre & Zillgitt, 1953; Metzger & Tubbs, 1971) culminated in development of a popular residual structure marking guide manual (Arbogast, 1957) which, along with other resources (e.g., Tubbs, 1977b), has since dominated the application of selection management (Kern et al., 2014). Under STS, select trees from all size classes are removed in harvests every 10 – 20 years (Neumann, 2015), removing unhealthy or ill- formed smaller diameter trees alongside larger, economically valuable stems. Single-tree gaps generate low-light understory environments (Beaudet et al., 2004) intended to promote sugar maple (Acer saccharum Marsh) and other shade-tolerant tree species (Crow et al., 2002; Angers et al., 2005; Poznanovic et al., 2013). Single-tree selection theoretically provides many benefits, including a steady stream of timber products and a low-impact harvesting system, which is visually appealing. However, recent studies suggest single tree selection may be yielding undesirable tree regeneration patterns, including lower diversity (Crow et al., 2002; Neuendorff et al., 2007; Powers & Nagel, 2009).and low density of desirable species, particularly sugar maple (Leak, 2006; Neuendorff et al., 2007; Powers & Nagel, 2009; Matonis et al., 2011). Furthermore, few 2 studies have analyzed age structure in selection-managed northern hardwood forests, instead relying on diameter distribution as a proxy; but age is vital to assess whether recruitment failures are recent or chronic. Given predicted declines in tree fitness (B. M. Rogers et al., 2017) and climatic changes (Byun & Hamlet, 2018) in the Great Lakes region, characterizing patterns of tree regeneration, recruitment, and associated potential driving factors is vital for promoting long-term forest health and sustainability via management. Management strategies are interacting with a variety of abiotic and biotic factors to influence forest dynamics. Forests in the Great Lakes are fundamentally distributed by post- glacial landforms and their associated nutrient availability and soil-water holding capacity (Zak et al., 1986, 1989; Baribault et al., 2010). Underlying soil type is further refined by a three-fold gradient of annual snow fall driven by the lake-effect, influencing tree species distribution and abundance (Henne et al., 2007). Northern hardwood forests vary in composition across the subset of soil types which support these forests. Interacting with the abiotic landscape, a variety of biotic drivers further filter species, including preferential browsing by white-tailed deer (Odocoileus virginianus) (Curtis & Rushmore, 1958; Horsley et al., 2003; Côté et al., 2004; Kain et al., 2011; Matonis et al., 2011; White, 2012; Bradshaw & Waller, 2016) and competition from woody non-tree shrubs (Royo and Carson, 2006; Walters et al., 2020a) and sedge (Randall & Walters, 2019). Invasive pests and pathogens have altered canopy composition and competitive outcomes for a variety of species, including American beech (Fagus grandifolia), white ash (Fraxinus americana), eastern hemlock (Tsuga canadensis), and American elm (Ulnus americana) (Parker & Leopold, 1983; Forrester et al., 2003; Nuckolls et al., 2009; Klooster et al., 2014). Although causal mechanisms between driving factors and northern hardwood forest dynamics are well researched, few studies have analyzed multiple drivers, across the range of 3 their values, and their associations with northern hardwood forest regeneration structure and outcomes. The goal of this dissertation is to characterize regeneration and recruitment outcomes in selection managed northern hardwood forests, with particular focus on a key species, sugar maple. Here, we use regeneration to refer to establishment of new tree seedlings and recruitment to characterize the process of small seedlings surviving and growing into taller size classes which have largely bypassed deer browsing and shrub competition, and therefore have a greater likelihood of reaching the overstory, compared to shorter classes. We assess regeneration and recruitment largely by stem densities. This dissertation relies on data from 141 selection- managed northern hardwood forests across northern Michigan, which cover a gradient of deer density, site quality, and canopy structure and composition. These stands are part of a larger research project analyzing regeneration outcomes to alternative silvicultural systems (Walters et al., 2020). This dissertation focuses on pre-harvest vegetation data and stand dynamics to characterize outcomes of decades of selection management in influencing forest structure and function. In Chapter 2 (Complex drivers of sugar maple (Acer saccharum) regeneration reveal challenges to long-term sustainability of managed northern hardwood forests), I address one of the most concerning trends to forest management: declining sugar maple regeneration patterns. Exploring this question through a Bayesian modeling framework, I analyzed how a variety of plot and stand level factors influenced sugar maple regeneration outcomes, drawing on vegetation surveys collected in all 141 stands. This chapter tests the hypotheses that sugar maple regeneration will: 1) be highly variable, with several regions of unacceptably low sapling densities; 2) negatively associate with historical deer use, particularly for sapling layers in the 4 deer browsing zone; 3) increase with stand-level density of mature, seed-bearing sugar maple trees, particularly for seedlings; 4) decrease with increasing woody non-tree vegetation density for the seedling class; 5) decline with higher total canopy density; and 6) increase with higher site quality. In Chapter 3 (Sugar maple age structure in northern hardwood stands managed by selection silviculture), I analyzed whether current low regeneration densities of sugar maple reflect a recent regeneration and recruitment failure or whether recent patterns reflect a chronic issue. To address this question, I analyze basal discs from recently harvested sugar maple stems for a subset of 51 stands, which spanned the extent of the larger study as well as landscape and stand-level drivers. With a collection of 1499 sugar maple tree samples, I addressed the following hypotheses: 1) Stands will be dominated by a canopy cohort of sugar maple stems, approximately 100 – 120 years old; 2) Stands will have one or two main age cohorts (less than the expected four age cohorts) which are > 60 years old; 3) Saplings (5 – 10 cm) will be older (> 60 years) than past studies indicate; 4) Older, sparser saplings will associate with areas of historic high deer density and/or lower site quality; and 5) Age and diameter will have a non- linear relationship, with a wide range of tree diameters having similar ages. In Chapter 4 (Characterizing seedling and sapling species diversity in managed northern hardwood forests portends low stand resiliency), I considered forest resiliency in a broader context, analyzing stand-level tree diversity with a focus on regeneration. In this analysis, I characterized how diversity varies with size class and with key landscape and stand-level drivers. To infer greater nuance to diversity measures, I also analyzed how individual tree species respond to the same suite of drivers for key regeneration size classes. This chapter addresses the following hypotheses: 1) Tree diversity will be highest in the canopy, followed by the seedling 5 class, and lowest for saplings class, and will be higher on lower quality versus higher quality sites; 2) Canopy composition will associate with regeneration density and diversity via seed source limitation, particularly for large-seeded species; 3 )Highly palatable tree species, like sugar maple, will negatively associate with deer use for the sapling class and will have greater representation in the canopy or seedling layers versus the sapling layer; 4) Shade intolerant and mid-tolerant species will associate negatively with increasing canopy basal area and will have greater representation in the canopy versus the understory; and 5) Small seeded species will face seedling establishment limitations, negatively associating with hardwood litter coverage, and therefore having greater representation in the overstory than in the understory. In Chapter 5 (Outcomes and implications of 60+ years of single-tree selection management in northern hardwood forests), I summarized the overarching findings of my dissertation. I explore my results in light of implications for single-tree selection management, climate change, and ongoing challenges to the long-term health and functioning of northern hardwood forests. 6 CHAPTER 2 COMPLEX DRIVERS OF SUGAR MAPLE (ACER SACCHARUM) REGENERATION REVEAL CHALLENGES TO LONG-TERM SUSTAINABILITY OF MANAGED NORTHERN HARDWOOD FORESTS 2.1 Abstract Single-tree selection silviculture management of northern hardwood forests relies on natural tree regeneration for long-term sustainability, yet current trends in tree regeneration and recruitment elicit concern. Low densities of economically valuable sugar maple (Acer saccharum Marsh.) in understories are often common, likely driven by many factors, including deer browsing, management-dictated stand structure, and site fertility/moisture regimes. However, landscape sugar maple regeneration patterns and relationships with underlying factors are largely unknown. We quantified associations of spatially varying factors with sugar maple regeneration using detailed vegetation and white-tailed deer winter fecal pellet surveys from 141 northern hardwood stands in Michigan, managed for decades with single-tree selection silviculture. We developed models of plot-level sugar maple regeneration counts for three key size classes as a function of plot- and stand-level predictors, including deer use, forest structure, and site quality. Among our 141 stands, sugar maple seedlings (< 50 cm tall) were consistently abundant, averaging 69,000 stems ha-1 and present in 76% of plots per stand, on average (25 plots per stand, each 2 m2). In contrast to seedlings, small (50 - 137 cm tall) and large (> 137 cm tall and < 5 cm DBH) sugar maple sapling densities were much lower, averaging 2,300 and 1,100 stems ha-1, respectively, and occurring, on average, in 31% and 32% of plots per stand (25 plots per stand, each 12.6 m2). Under a wide range of potential sugar maple stocking criteria, most stands are understocked in saplings, e.g., only 29% of stands had > 2,500 small sugar maple sapling stems ha-1 and only 7 35% of stands had > 1,000 large sapling stems ha-1. Based on our models, sapling densities negatively associated with deer use and were most abundant on medium quality sites. Across all size classes, negative associations with subcanopy trees and/or shrub densities suggest light limitation, whereas positive associations with sugar maple canopy trees > 25 cm DBH suggest persistent seed limitations. Overall, our study supports a need for alternative forest and/or deer management strategies over much of the range of northern hardwood forests in Michigan to promote higher densities of large sugar maple regeneration for canopy recruitment. Medium quality sites with abundant large sugar maple canopy trees and low deer browse pressure (for example, the deep snow region in the northern portion of the western Upper Peninsula) are the exception; high sugar maple sapling densities suggest forests in this region are thriving under single-tree selection management. 2.2 Introduction 2.2.1 Management and regeneration dynamics Sustainable forestry depends on regeneration and canopy recruitment of tree species favored by management. However, temperate forests globally exhibit recent trends of declining regeneration (Ramirez et al., 2018; Miller & McGill, 2019). In the northern United States, many naturally regenerated forests are understocked in regeneration size classes, and most are projected to undergo species compositional shifts (Vickers et al., 2019). Northern hardwood forests, which occur from the Great Lakes region to Maritime provinces (Braun, 1950), are one such naturally regenerated forest type with regeneration and recruitment concerns (Jenkins, 1997; R. O. Miller, 2004; Donovan, 2005). Seedling and sapling layers in managed northern hardwood forests may have lower species diversity compared to unmanaged stands (Neuendorff et al., 2007; Powers & Nagel, 2009), which undermines long-term forest resilience (Thompson et 8 al., 2009). Declining density of sugar maple (Acer saccharum Marsh.) regeneration throughout much of its northern range (Leak, 2006; Neuendorff et al., 2007; Powers & Nagel, 2009; Matonis et al., 2011) is particularly concerning as it currently dominates the canopy (Schulte et al., 2007), is important ecologically as a late successional species (Curtis & McIntosh, 1951), and is highly valued economically (Linehan & Jacobson, 2005; Duval et al., 2014). Given that northern hardwood management systems are designed to promote abundant sugar maple regeneration and recruitment, understanding stand and landscape level factors predictive of sugar maple regeneration that managers can manipulate is key for promoting long-term sustainability. In the Great Lakes region, low sugar maple regeneration densities exist following decades of single-tree selection management. Many of these second growth stands originated 100+ years ago as primarily even aged, following widespread and intense logging in the 19 th and early 20th centuries (Dickmann & Leefers, 2016). Since stand origination, species richness and structural diversity have declined, in part driven by natural succession and by selective removal of species for economic reasons (Nyland, 1992; Schulte et al., 2007). Single-tree selection was introduced in the 1950’s to convert a preponderance of relatively young, even-aged northern hardwood forests to uneven aged stands and/or to maintain uneven-age structure over the limited extents this structure occurred (Arbogast, 1957). In theory, this regime involves partial stand harvesting of dispersed, individual trees from all size classes, generating small gaps at 10–20- year harvest intervals (Neumann, 2015). In practice, single-tree selection prevails within stands but with an increased emphasis in the last couple of decades on including some small group harvest areas within stands, which has been implemented to varying degrees. Given the dominance of single-tree selection within and across stands, we hereafter use the term single-tree selection. Regardless of specific pattern, stand basal area (BA) is typically reduced from ~ 25 - 9 30 m2 ha-1 to 17 - 21 m2 ha-1 (Pond et al., 2014; Neumann, 2015; OMNRF, 2015). In theory, this mimics frequent, low-intensity disturbances characteristic of unmanaged northern hardwood stands (Frelich & Lorimer, 1991) but contrasts with high intensity disturbances that led to current forest origination (Whitney, 1987). With consistent regeneration and recruitment, single-tree selection generates a continuous stream of timber products, offering economic stability. Since introduction, single-tree selection has been the dominant management paradigm for northern hardwood forests, particularly in the Great Lakes region (Kern et al., 2014). Lower-light conditions characterizing single-tree selection (e.g., at 5 m above the ground, ~10-20% of above-canopy photosynthetic photon flux density for the first few years following harvest, (Beaudet et al., 2004) are intended to promote sugar maple, a shade-tolerant species (Tubbs, 1977a), but recent findings challenge this notion. Harvest regimes which generate low (e.g., single-tree selection) and high (e.g., clearcutting) light conditions tend to have lower sugar maple regeneration densities compared to harvests of intermediate canopy openings (e.g., 50% basal area reduction) and light availability (e.g., ~80% of ambient photosynthetically active radiation for shelterwood cutting, Grayson et al., 2012) (Matonis et al., 2011; Cleavitt et al., 2018; Danyagri et al., 2019). Instead of creating long-term light conditions suitable for sapling recruitment, light levels in single-tree canopy gaps may decline via horizontal expansion of gap edge canopy trees faster than advanced regeneration can recruit (Kern et al., 2013), with light levels returning to pre-harvest conditions within 8-11 years (Beaudet et al., 2004). The notion that small gaps are insufficient for recruitment is supported by observation that successful sugar maple recruitment was more likely in larger single-tree gaps (~78 m2) with taller (> 4 m) advanced regeneration (Cole & Lorimer, 2005). Since sugar maple saplings decline in shade tolerance and growth rate as their size/age increases (Donoso et al., 2000; Sendall et al., 2015), 10 vertical growth may stall and mortality increase unless light levels are maintained or increased in subsequent single-tree selection harvests. Over several cutting cycles, use of harvest gaps too small to recruit saplings to the overstory may lead to a surplus of old, suppressed sugar maple and other shade-tolerant species accumulating in the understory. With limited recruitment from sapling classes over several decades of continued single-tree selection management, canopy/subcanopy tree classes become understocked, first in the smaller classes and ultimately in the largest classes (Millington et al., 2011). Consistent with limited recruitment over decades, but not demonstrating cause and effect, low small tree (pole class 10 – 25 cm DBH) density in northern hardwood stands has been identified in single-tree selection managed northern hardwood stands in northern Michigan (Walters et al., 2020b) and broadly at the landscape level (Hanberry & Abrams, 2019). 2.2.2 Local and landscape factors affecting regeneration Seedling and sapling competition for light is not only influenced by distant overstory canopies but also by understory shrubs and larger sapling classes (Schwinning & Weiner, 1998; Collin et al., 2017). Shrubs broadly impact regeneration outcomes and shift forest stand dynamics (Royo & Carson, 2006). For sugar maple, shrubs inhibit seedling growth, particularly on high quality sites, though they can improve survival during periods of drought (Berkowitz et al., 1995). Dense, tall tree regeneration can also cast deep shade on smaller regeneration below, inhibiting growth and survival even for shade-tolerant species such as sugar maple (Beaudet et al., 2002). For example, American beech (Fagus grandifolia Ehrh.) can establish dense sapling layers that associate with low densities of subordinate sugar maple stems, suggesting competitive effects (Hane, 2003; Nyland et al., 2019; Elenitsky et al., 2020); beech's competitive edge is likely exacerbated by beech bark disease, which triggers dense beech sapling thickets (Forrester 11 et al., 2003). For single-tree selection managed forests, shrubs and taller, dense regeneration layers present potential competitive challenges for sugar maple regeneration. Fundamentally, sugar maple seedling regeneration is driven by seed availability and site quality, which affects germination, growth, and survival. Seedling density is closely tied to seed production (Bjorkbom, 1979), with larger trees producing more seed (Garrett & Graber, 1995). Sugar maple is also known to regenerate vigorously via stump sprouting after harvest (Forrester et al., 2014). Sugar maple trees (and therefore seed source/stump sprouting potential) are generally found on high nutrient and moderate to high water holding capacity soils (Burger & Kotar, 2003). Key resources potentially varying and limiting the establishment, growth, and survival of sugar maple regeneration include water, nitrogen, and calcium. Variation in some or all of these resources likely underlie variation in site quality and have been found to scale with ecological classification systems (Zak et al., 1989) and habitat classification categories (Walters & Reich, 1997; Baribault et al., 2010). High quality site conditions promote germination success (Tubbs, 1977a), enhance seedling growth (Walters & Reich, 1997), and reduce regeneration mortality rates (Burns & Honkala, 1990; Caspersen & Kobe, 2001; Kobe, 2006), meaning sugar maple has previously and is expected to thrive on high quality sites. Northern hardwood regeneration exhibits differential abundance among site qualities (Elenitsky et al., 2020), though associations can be confounded at smaller spatial scales (Matonis et al., 2011). Since sugar maple canopy growth rates are higher on higher quality sites (Baral et al., 2016) and individual super-producing masting trees tend to associate with higher nutrient availability (Minor & Kobe, 2017), high quality sites may also have greater seed production, leading to compounding effects of site quality and seed source on seedling density. 12 Additionally, abundant populations of white-tailed deer (Odocoileus virginianus Zimmermann, hereafter called deer) influence regeneration patterns. At broad scales, deer are linked to understocked regeneration throughout northern United States (Vickers et al., 2019) and shifting tree regeneration and herbaceous communities via preferential browsing (Horsley et al., 2003; Côté et al., 2004; Bradshaw & Waller, 2016). Deer exclusion experiments in northern hardwood forests have quantified regeneration density and average height reductions caused by deer (Curtis & Rushmore, 1958; Horsley et al., 2003; Kain et al., 2011; White, 2012). Regarding tree regeneration, high deer densities can mitigate positive effects of disturbance (Nuttle et al., 2013), increase optimal gap size (Walters et al., 2016), and shift shrub coverage from competitive to facilitative (Walters et al., 2016). Browsing pressure particularly suppresses small saplings in the browse at-risk zone (<2 m in height; (Walters et al., 2020a), inhibiting recruitment. While direct impacts of deer on vegetation are well-documented, relative importance of continuous variation in browsing pressure over regional landscapes, with other factors also varying, is largely unknown. Quantifying landscape level patterns of sugar maple regeneration and identifying predictive factors for managed northern hardwood forests is imperative for continued sustainable management; long-term regeneration failure has enormous impacts, economically and ecologically. Although likely mechanisms from small-scale studies are useful, these studies lack predictive power across broad areas. In contrast, geographically widespread studies using public datasets help identify broad patterns but lack detail and rigor to assess all potential predictive factors. In this study, we bridge the gap between geographic breadth and rigorous detail to characterize landscape patterns and identify factors associated with sugar maple regeneration in single-tree selection managed northern hardwood forests, across northern Michigan. We expect 13 that sugar maple regeneration will: 1) be highly variable, with several regions of unacceptably low sapling densities; 2) negatively associate with deer use, particularly for sapling layers in the deer browsing zone; 3) increase with stand-level density of mature, seed-bearing sugar maple trees, particularly for seedlings; 4) decrease with increasing woody non-tree vegetation density for the seedling class; 5) decline with larger total canopy density; and 6) increase with higher site quality. We interpret analysis results with consideration of management implications and highlight conditions and areas that are particularly favorable or unfavorable for continuing use of single-tree selection. 2.3 Methods 2.3.1 Study area We quantified vegetation and deer use in 141 managed northern hardwood stands throughout northern Michigan (Walters et al., 2020b). At the start of the study, stands were considered ready for partial harvest by standard single-tree selection criteria (i.e., > 23 m2 ha-1 BA and well stocked in sawtimber classes; Arbogast, 1957). Most of the 141 stands were State owned (n=119) and managed by the Michigan Department of Natural Resources (MDNR), with private industrial sites (Hancock Timber Resource Group, The Rohatyn Group, n=22) comprising the remaining stands. Stands were generally on upland sites and dominated by sugar or red maple (Acer rubrum L.). In addition to sugar and red maple, other common northern hardwood tree species included striped maple (Acer pensylvanicum L.), balsam fir (Abies balsamea [L.] Mill.), yellow birch (Betula alleghaniensis Britton), paper birch (Betula papyrifera Marsh.), American beech, white ash (Fraxinus americana L.), ironwood (Ostrya virginiana Mill.), white spruce (Picea glauca [Moench] Voss), white pine (Pinus stobus L.), bigtooth aspen (Populus grandidentata Michx.), white oak (Quercus alba L.), red oak (Quercus 14 rubra L.), basswood (Tilia americana L.), hemlock (Tsuga canadensis L.), and American elm (Ulmus americana L.). Northern hardwood forests of Michigan are broadly found on mesic to wet mesic sites with fertile medium-textured upland soils (Dickmann & Leefers, 2016). In the Great Lakes region, nutrient availability and soil moisture tend to positively covary (Zak et al., 1986; Walters & Reich, 1997), and adjacent Great Lakes drive significant patterns in snowfall and overall annual precipitation. Together, these factors are strong drivers of northern hardwood forest distribution (Burger & Kotar, 2003; Henne et al., 2007). Our study design was intended to maximize variation in site quality among stands (Figure 2.1), which we characterized with a habitat classification system (Burger & Kotar, 2003). 15 Figure 2.1. Locations of 141 single-tree selection-managed northern hardwood stands sampled for sugar maple regeneration and recruitment across northern Michigan by site class. SQ1 stands are poor to poor/medium, SQ2 medium, SQ3 medium/rich, and SQ4 rich to very rich (Burger and Kotar, 2003). 2.3.2 Field methods In summer 2016, we established a permanent 12.14 ha square within each of our 141 sites, which we heretofore refer to as stands, with 25 survey points laid out in a systematic grid; survey points were spaced 69.5 m apart in the grid. We surveyed vegetation in these stands in summer 2017. At each survey point, we established two 1 m2 quadrats centered 3 m east and west of the survey point, which together represent one plot (2 m2). In those quadrats, we counted all tree seedlings, by species, shorter than 50 cm, and ocularly estimated percent cover of seedlings and shrubs. Using the survey point as a common center, we used a 2 m radius circular plot (12.6 m2) to tally tree regeneration (by species) 50 cm to 137 cm tall (small saplings), and tally and measure DBH (by species) for stems > 137 cm tall and < 5 cm DBH (large saplings). In a 6 m radius circular plot (113 m2), we tallied and measured DBH for stems > 5 cm DBH. As we 16 traversed stands, we recorded herbaceous plant assemblages and key indicator species to assign a dominant site-quality index (Burger & Kotar, 2003) at the stand level. To assess winter deer use, we conducted pellet surveys in spring of 2017 at a subset of 50 sites, stratified spatially to sample our entire study region. In spring of 2019, entering the first growing season following timber harvest (unrelated to this analysis) of the sites, we surveyed 139 sites using the same protocol. Although pellet surveys are error-prone, they are cost- effective, feasible at large scales, and provide a reasonable approximation of winter deer use (Forsyth et al., 2007; Urbanek et al., 2012; Brinkman et al., 2013). We surveyed linear transects (6 m wide) that totaled 628 m in length, spatially dispersed throughout the stand (Appendix A). Transects were subdivided into 103 segments ~ 6 m long in which we recorded deer pellet occupancy. We conducted surveys from 22 April to 20 May in 2017 and 14 April to 29 May in 2019, prior to emergence of spring ephemerals, and assumed leaf-off to be 01 November for all stands, for both winters. With pellets atop last-years leaf litter, the method assumes pellets were deposited between leaf-off and date of survey. We developed a stand-level deer use model which incorporated climatic and landcover predictors to provide point estimates of deer use for stands which were not surveyed in 2017 (methods in Appendix A). 2.3.3 Sugar maple regeneration density modeling framework Our goal was to estimate parameter values and associated uncertainty for a set of carefully selected ecological factors potentially predictive of sugar maple regeneration count data in key size classes. All analyses and modeling were conducted in R (R Core Team, 2018). We developed a multilevel model (predictors at plot- and stand-level, Table 2.1) to predict plot-level sugar maple counts, which we ran separately for sugar maple counts in three size classes: seedlings (0 – 50 cm tall), small saplings (50 cm – 137 cm tall), and large saplings (> 137cm tall 17 and less than 5 cm DBH). These size classes represent key stages of sugar maple size and canopy recruitment probability with respect to deer browsing pressure (active driver for small saplings, potential legacy effects on large saplings escaped from the deer browse zone, and less tied to seedlings), seed source (anticipated greatest association with seedlings). Large saplings are considered to have highest potential for canopy recruitment, seedlings the lowest. Plot sizes were 2 m2 for seedlings and 12.6 m2 for small and large saplings. Plot-level predictors in our model included percent shrub coverage (SHRUB, all non-tree woody species, generally < 100 cm tall in our stands) and total stem count of relatively taller regeneration classes, all species (SAP1, small saplings; SAP2, large saplings; and SUB, subcanopy trees, 5 – 10 cm DBH). Stand level predictors included modeled average estimated likelihood of deer use (DEER, see Appendix A for deer modeling details; Figure 2.2), four ordinal categories of site quality (SQ1-4, 1 being poorest, see Appendix B for details), average total BA trees > 10 cm DBH (BA, proxy for forest canopy effect on light availability and microclimate), and average BA of sugar maple trees > 25 cm DBH (SMBA, proxy for seed availability, size criteria determined using methods in Appendix C). BA and SMBA were averaged at the stand level to better reflect management scale and because a 6m radius plot may not fully capture canopy competition or seed availability. 18 Name Description Mean (Range) Plot SAP1 Total stem counts of small saplings (> 50 cm tall 6 (0 – 116) level and < 137 cm tall) per plot (12.6 m2), all species SAP2 Total count of large saplings (> 137 cm tall and < 5 4 (0 – 38) cm DBH) per plot (12.6 m2), all species SUB Total count of subcanopy trees (5 – 10 cm DBH) 3 (0 – 33) per plot (12.6 m2), all species SHRUB Ocular estimate of percent ground obstructed by 2 (0 – 93) shrubs, viewed from above (0 – 100%) Stand BA Basal area all species > 10 cm DBH (m2ha-1) 26 (19 – 44) level SM.BA Basal area of sugar maple > 25 cm DBH (m2ha-1) 13 (0 – 27) DEER Modeled deer use (0 – 100%) 25 (4 – 76) SQ1 Poor to poor/medium site quality (n = 39) - SQ2 Medium site quality (n = 22) - SQ3 Medium/rich site quality (n = 47) - SQ4 Rich to very rich site quality (n = 33) - Table 2.1. For our multilevel models of the plot-level count of sugar maple regeneration (seedlings, small saplings, and large saplings), we included predictor variables at both plot and stand levels. The table lists a brief description, the abbreviation used, the mean, and the range of parameter values. 19 Figure 2.2. Model estimated stand-level deer use, winter of 2017. Zero represents a stand with no estimated winter deer use (0 segments had fecal pellets), while 100 indicates a high use stand where all 103 transect segments (each segment ~6m x 6m) had fecal pellets. We checked for interactive effects between site quality and three broad regions of our study (Elenitsky et al., 2020) and between site quality and deer use (selective browsing pressure and/or deer occupancy within region or site class may be altered via differences in species assemblages of potential browse species); analyses did not support adding any interactive effects (Appendix C). Multicollinearity among the predictor variables was assessed using variance inflation factors (via vif in the car package (Fox & Weisberg, 2011)); at a threshold of 2, none warranted removal. To understand how stands covary in their characteristics (e.g., do high quality stands tend to have more canopy sugar maple trees?), we calculated a correlation matrix based on Kendall’s tau statistic (n=141), averaging plot-level predictors to analyze all at the stand level; p-values (< 0.05 significance threshold) were calculated using cor_pmat from the rstatix package (Kassambara, 2020). We also calculated summary statistics to characterize our 20 response data (sugar maple regeneration), including occurrences, average densities, relative abundance, and boxplots of regeneration density vs. site quality. For each regeneration size class, we developed five candidate models (Appendix C for details on model formulation and associated diagnostic checks), including a null model and various inclusions of quadratic terms to account for non-linearity, a common approach in regeneration modeling (e.g., Horsley et al., 2003; Schwarz et al., 2003; Danyagri et al., 2019; specifically for deer use effects: Horsley et al., 2003; Tremblay et al., 2006; Bradshaw & Waller, 2016). We selected models based on deviance information criterion (DIC) comparison, which reflects model complexity and fit (Spiegelhalter et al., 2002). For each size class, we used the following model specification to estimate sugar maple stem counts y in stand i and plot j: Plot level 𝑦𝑖𝑗 ~ 𝑁𝑒𝑔𝐵𝑖𝑛(𝛼𝑖 + 𝛽𝑝 𝑿𝒑𝒊𝒋 , δ𝑖 ) Stand level 𝛼𝑖 ~ 𝑁𝑜𝑟𝑚(𝜇 + 𝛽𝑠 𝑿𝒔𝒊 , 𝜎 2 ) 𝑦𝑖𝑗 is sampled from a negative binomial distribution with mean 𝛼𝑖 + 𝛽𝑝 𝑿𝒑𝒊𝒋, where 𝛼𝑖 is a stand-specific intercept and 𝛽𝑝 and 𝑿𝒑𝒊𝒋 represent the vector and model matrix for plot-level predictors SHRUB, SAP1 (for seedling model), SAP2 (for seedling and small sapling model), and SUB plus any included quadratic terms (Table 2.1). δ𝑖 is a stand-specific, plot-level dispersion parameter. 𝛼𝑖 is drawn from the normal distribution with mean 𝜇 + 𝛽𝑠 𝑿𝒔𝒊 , where 𝜇 is the grand mean and 𝛽𝑠 and 𝑿𝒔𝒊 represent the vector and model matrix for stand-level predictors BA, SM.BA, DEER, SQ1, SQ3, and SQ4 plus any included quadratic terms (SQ2 included in the intercept, Table 2.1). 𝜎 2 represents stand-level variation. We estimated stand-specific dispersion parameters to account for variation among stands in plot-level variance. We mean centered all predictors and divided by two standard deviations to enable comparison (Gelman & Hill, 2007). 21 Predictive inference was based on 10,000 post burn-in samples from three Markov chain Monte Carlo (MCMC) chains using the RJAGS package (Plummer, 2018), which is an interface to JAGS (Just Another Gibbs Sampler) software. We assessed model convergence using Gelman Rubin diagnostics (coda package), trace plots, and residuals vs. predicted plots. Histograms of predicted zeros compared to true data zeros for each model indicated that the models adequately predicted absences, indicating that the negative binomial was an appropriate fit. We checked for spatial autocorrelation in model residuals, both plots within individual stands and among stands, using variograms and found no residual spatial autocorrelation in the model (Appendix C). To visualize the model results, we simulated model predictions for 1,000 iterations across the range of predictor values for each significant continuous predictor, setting all other predictors at their mean and using average dispersion parameters and stand-level intercepts. To plot the model predictions, we applied a loess smoothing curve on predicted means and 95% lower and upper confidence intervals. 2.4 Results 2.4.1 Sugar maple regeneration characteristics Stand-level sugar maple regeneration was highly variable within and among size classes (Table 2.2). For all three size classes, frequency distributions of plots with sugar maple were right skewed with wide ranges (Table 2.2), indicating that most plots had few or no stems and a small percentage had high stem counts (Table 2.2). Sugar maple seedlings (0 – 50 cm tall) were widespread (present in 99% of stands and 76% of all plots) and abundant (average 14 stems per 2 m2 plot, or 70,000 stems ha-1) (Table 2.2). Small and large sugar maple saplings were also widespread, occurring in 79 and 90% of stands but only occurred in 31 and 32% of total plots, respectively (Table 2.2). In 8% of our stands, our surveyed area failed to capture either a single 22 small or large sugar maple sapling. Average counts for small and large saplings were 3 and 1 stems per 12.6 m2 plot, respectively, or approximately 2300 and 1090 stems ha-1. Median densities for small and large sapling stems were lower than the mean at 600 stems ha-1 and 500 stems ha-1, respectively (Table 2.2). We observed regional variation in sugar maple stand sapling densities (seedling densities were more uniformly dispersed, Figure 2.3A). For small saplings, the northern Lower and southern half of Upper Peninsula had low densities, with values ≥ 10,000 stems ha-1 only found in the northern half of the Upper Peninsula (Figure 2.3B). Large sapling followed a similar, but somewhat weaker pattern as higher densities did exist in some stands in the southern half of the Upper Peninsula as well as in the Lower Peninsula (Figure 2.3C). 23 Plot level Stems per plot % relative abundance Plot occurrences Mean Median Range Mean Median Range Seedlings 14 5 0 – 245 60 75 0 - 100 2675 / 3524 (76%) Small 3 0 0 – 100 32 0 0 - 100 1082 / 3524 (31%) saplings Large 1 0 0 - 34 30 0 0 - 100 1130 / 3524 (32%) saplings Stand Thousand stems per ha % relative abundance Stand occurrences level Mean Median Range Mean Median Range Seedlings 69 47 0 – 434 64 71 0 - 99 140 / 141 (99%) Small 2.3 0.6 0 – 22 34 25 0 - 98 111 / 141 (79%) saplings Large 1.1 0.5 0–7 31 20 0 - 100 127 / 141 (90%) saplings Table 2.2. The mean, median, and range of sugar maple regeneration densities for seedlings, small saplings, and large saplings at the plot and stand level. Seedlings are 0-50 cm tall, small saplings are 50 cm - 137 cm, and large saplings are 137 cm tall or greater and up to 5 cm DBH. Densities are presented in stems per plot (2 m2 plot for seedlings, 12.5 m2 plot for small and large saplings). Plot occurrences indicates the number of surveyed plots which contained at least one sugar maple. Stand occurrences indicates the number of sites which contained at least one sugar maple stem in any of the plots surveyed in that stand. The percent relative abundance refers to the percentage of sugar maple regeneration stems relative to the total density of stems in that size class; for plot percentage relative abundance, plots with zero total stems in a size class were omitted. 24 A B C Figure 2.3. Average stand-level density of sugar maple seedlings (0 – 50 cm tall; A), small saplings (50 – 137 cm tall; B), and large saplings (>137 cm tall and <5 cm diameter breast height) based on data collected for 141 northern hardwood stands, summer 2017. 25 2.4.2 Stand characteristic covariation Among variables used to predict sugar maple regeneration densities, we found that deer use, shrub cover, sugar maple canopy BA, and site quality positively covaried in our stands (Figure 2.4). These variables generally negatively covaried with density of total understory vegetation (including large saplings and subcanopy trees). Additionally, total stand BA positively covaried with canopy sugar maple BA and negatively covaried with shrub cover. Figure 2.4. Kendall’s tau correlation matrix of predictor variables used to model stand level sugar maple regeneration for Michigan (non-significant (p>0.05) correlations are labeled with an x, positive associations are blue, and negative associations are red with parallel lines). Variables include total stem density of small saplings (SAP1), total density of large saplings (SAP2), total density of trees 5 – 10 cm DBH (SUB), deer usage (DEER), site quality (SQ, treated as continuous from 1-4), shrub coverage (SHRUB), basal area of sugar maple trees > 25 cm (SM.BA), and total stand basal area (BA). 2.4.3 Sugar maple regeneration models For sugar maple seedling counts, the model intercept was ~2.5 on the log scale, predicting occurrence more often than absence. At the plot level, seedling counts negatively associated with shrub cover, and exhibited quadratic relationships with trees 5- 10 cm DBH, and 26 small and large sugar maple saplings (Figure 2.5). Sugar maple seedlings generally increased with total small sapling counts (eventually plateauing) but generally decreased with total large saplings counts, total sub-canopy tree counts, and percent shrub cover (Figure 2.6). At the stand level, sugar maple seedling count increased with BA of canopy sugar maple trees until ~18 m2 ha-1, then declined (Figure 2.6). There were no clear overarching associations with site quality (Figure 9), though stand-level densities of sugar maple seedlings on medium/rich sites (SQ3) were significantly lower than on medium sites (SQ2, included in the intercept) (Figure 2.5). Deer use and total stand basal area were not significantly associated with seedling counts. Our top- ranked sugar maple seedling model ranked higher than the null model (Appendix C). 27 Figure 2.5. Model parameter estimates for predictors of sugar maple regeneration, with 95% Bayesian credible intervals. Variables include total stem density of small sugar maple saplings (SAP1), total density of large sugar maple saplings (SAP2), total density of trees 5 – 10 cm DBH (SUB), shrub coverage (SHRUB), deer usage (DEER), total stand basal area (BA), basal area of sugar maple trees > 25 cm (SM.BA), site quality (SQ1 is poor, SQ3 medium/rich, SQ4 rich to very rich, and SQ2 medium included in the intercept), and model intercept (MU), which represents average plot-level stem count, on the log scale, when all predictors are at average values. Quadratic terms (indicated by superscripted 2) included in the final models portrayed with dashed lines. 28 Figure 2.6. Model estimates for the relationship of sugar maple seedlings per plot (2 m2) with the plot-level densities of small saplings (SAP1), large saplings, (SAP2) and all sub-canopy trees 5 – 10 cm DBH (SUB); the stand-level basal area of sugar maple canopy trees larger than 25 cm DBH (SM.BA); and the plot-level density of shrubs (SHRUB). Mean (black line) and 95% Bayesian credible intervals (gray lines) are displayed, based on the average dispersion parameter estimate. Small and large sugar maple sapling counts associated similarly with predictor variables. Overall intercepts were not significantly different from zero, meaning both models more often predicted absence than a count (Figure 2.5). Counts of small and large sugar maple saplings decreased with increasing deer use, increased with basal area of sugar maple trees > 25 cm DBH at lower basal areas, and were unassociated with overall stand BA (Figures 2.5, 2.7, 2.8). Sugar maple saplings were most abundant on low (SQ1) or medium quality sites (SQ2, included in the intercept) and lowest on medium/rich (SQ3) and rich sites (SQ4) (Figure 2.9). Small sugar maple saplings were positively related to total large saplings (weakly at low density), and large sugar maple saplings showed no association to subcanopy trees (Figure 2.7). Similar to seedlings, 29 small saplings declined with increasing counts of subcanopy trees. Small and large sugar maple saplings differed in that shrub density was weakly negatively related to large sapling counts (Figure 2.8) and unassociated with small saplings. Both models ranked higher than the null model (Appendix C). Figure 2.7. Model estimates for the relationship of sugar maple small saplings per plot (12.6 m2) with large saplings per plot (SAP2), subcanopy trees per plot (SUB), deer use (DEER), and basal area of sugar maple canopy trees larger than 25cm DBH (SM.BA). Mean (black line) and 95% Bayesian credible intervals (gray lines) are displayed for the average dispersion parameter estimate. 30 Figure 2.8. Model estimates for the relationship of sugar maple large saplings per plot (12.6 m2) with deer use (DEER), basal area of sugar maple canopy trees larger than 25cm DBH (SM.BA), and shrub coverage (SHRUB). Mean (black line) and 95% Bayesian credible intervals (gray lines) are displayed for the average dispersion parameter estimate. 31 Figure 2.9. For each sugar maple regeneration class, boxplots of stand-level average density by site quality category (sensu Kotar and Burger 2003), northern Michigan, 2017. Site quality categories include poor (SQ1), medium (SQ2), medium to rich (SQ3), and rich to very rich (SQ4) in the model. Middle lines represent median values, the lower and upper hinge represent the first and third quartiles, and the whiskers represent values no further than 1.5 of the interquartile range (distance from first to third quartile); individual points are outliers. 32 2.5 Discussion Widespread and growing evidence suggests understocked sugar maple regeneration in single-tree selection-managed northern hardwood forests (Jenkins, 1997; Miller, 2004; Donovan, 2005; Leak, 2006; Neuendorff et al., 2007; Powers & Nagel, 2009; Matonis et al., 2011), with uncertainties regarding full extent and associated predictive factors; this information is crucial for long-term management success. Detailed studies which have elucidated mechanism tend to be limited in geographic scope and stand characteristics, while broad-scaled analyses gain power but may lack the detailed information needed to identify ecological drivers. Drawing on established theory and previous insight, we bridged depth and breadth by surveying 141 single- tree selection managed northern hardwood stands across northern Michigan. Our results confirm widespread and concerning trends in sugar maple regeneration and recruitment. Sugar maple regeneration densities were variable for key size classes and low across wide geographic ranges and stand characteristics. While densities of sugar maple seedlings were generally high (though potentially limited on some sites by low mature sugar maple seed sources and/or competing taller vegetation), densities of small and large sugar maple saplings were overall low, particularly in portions of the northern Lower Peninsula and the central southern Upper Peninsula. Densities in the sapling classes were approximately three orders of magnitude less than for the seedling class. In 20% of our stands, our surveyed area (totaling 315 m2 per stand) did not include a single sugar maple small sapling (50 cm – 137 cm tall). We found evidence of potential sugar maple regeneration limitation by both deer browsing and seed availability in our stands. Negative associations of deer use with small saplings suggest proximal/recent deer browsing impacts; effects for large saplings that have largely escaped deer browse pressure may indicate long-term legacies of deer browsing 33 inhibiting sapling canopy recruitment, which aligns with previous studies (Horsley et al., 2003; Kain et al., 2011; Matonis et al., 2011; Bradshaw and Waller, 2016). We found no association between deer use and sugar maple seedling count, unsurprising given that shorter regeneration is often less browsed (perhaps protected by winter snow) and/or maintained in that size class by persistent browse pressure (Saunders & Puettmann, 1999; Randall & Walters, 2011). Our results also suggest persistent legacies of seed limitation in larger size classes. As expected (Bjorkbom, 1979), seedling count increased with sugar maple canopy BA (up to ~ 20 m2 ha-1); the decline of seedlings at higher sugar maple canopy BAs may point to potential light limitation at high densities of canopy sugar maple, as sugar maple trees cast deep shade (Canham & Burbank, 1994). However, general positive associations of small and large sugar maple saplings with sugar maple canopy BA suggest that despite other pressures, such as deer, stands with greater seed source tend to have more sugar maple saplings. Promoting abundance of sugar maple trees > 25 cm DBH, up to densities of 20 m2 ha-1, may lead to greater densities of sugar maple saplings even in regions of high deer populations. Light limitation from shrubs (Matonis et al., 2011; Kern et al., 2013) and larger tree regeneration (Schwinning & Weiner, 1998) influences regeneration outcomes in northern hardwood stands. Our results support previous work, with negative associations between sugar maple regeneration with competing larger regeneration (excluding the relative next largest size class) and/or shrub density. We attribute positive or neutral associations between a given sugar maple size class and total counts of the relative next largest size class (e.g., sugar maple seedlings positively associated with total small saplings, but negatively with large saplings and subcanopy trees) to categorizing a continuous distribution of regeneration heights for a species of generally high relative importance (also found in Elenitsky et al. 2020). We did not find evidence 34 supporting light limitation from canopy BA, previously identified as limiting to sugar maple regeneration (Matonis et al., 2011; Cleavitt et al., 2018; Danyagri et al., 2019). However, though not significant, all three mean stand BA parameter estimates were negative in our model; we may have failed to detect a significant association due to the distribution of stand BA included in our study with fewer stands at very high BA (range 19 – 44 m2 ha-1, mean 26 m2 ha-1). Stand BA also does not represent light reaching the forest floor well when subcanopy or understory vegetation is dense or when composition of canopies (with different canopy light transmissivities among species) vary among sites (Beaudet et al., 2004). In our study, sugar maple sapling densities diverged from traditional affinity for high quality sites (Burger & Kotar, 2003) (Figure 2.9). Although high quality sites generally enhance regeneration growth (Walters & Reich, 1997) and reduce mortality rates (Burns & Honkala, 1990; Caspersen & Kobe, 2001; Kobe, 2006), small and large sugar maple sapling abundance was generally lowest on rich sites. Seedlings, in contrast, did not demonstrate much significant variation by site quality; thus, abundant seedlings on high quality sites may not translate to abundant saplings. One possible explanation is higher quality sites had higher deer use (Figure 2.4), likely due to the spatial distribution of site qualities (e.g., in the western Upper Peninsula, site quality and deer use both go from low to high along a north-south gradient) (Figure 2.1). Deer and site quality were also confounding in Matonis et al. (2011) which comprised a subset of our study area (also in the western Upper Peninsula). Lower sapling counts on high quality sites could also be associated with deeper shade, driven by the higher leaf area index (LAI) generally associated with higher quality sites that is largely driven by species with high LAI on those sites (e.g., sugar maple) (Canham & Burbank, 1994; Fassnacht & Gower, 1997). More rapid sugar maple growth on high quality sites (Baral et al., 2016) could result in more rapid closure of small 35 canopy gaps, reducing canopy recruitment potential and increasing understory mortality. It is also possible that out categorical site quality indicators were too coarse to adequately capture variation in resources actually underlying seedling establishment, growth, and survival over sites (e.g., nitrogen, calcium, water); furthermore, site-quality as a stand-level predictor may fail to capture finer scale nutrient variability within stands which could affect sugar maple regeneration. While our study suggests that some combination of deer use, seed availability, light competition, and site quality drive sugar maple regeneration dynamics, our ability to draw conclusions on causal mechanisms is limited because our data are from a natural, snapshot in time experiment (Diamond, 1986). We addressed this potential pitfall by carefully choosing model predictors grounded in ecological theory. Due to the large number of stands and plots, our plot sizes were also relatively small, which may lower precision of our density estimates for a single plot. We overcame deer use sampling limitations by developing a deer use model (Appendix A) and using evidence-based statistical modeling approaches to account for non- linearities. Despite potential limitations, our study bridges the gap between broad-scale and detail-oriented, with a wide geographic extent and broadly varying site characteristics complemented by a rigorous and fine-scaled sampling design. This affords rare insight based on statistical power in understanding dynamics of managed northern hardwood forests. 2.6 Implications for Management Although Arbogast’s (1957) manual is considered definitive for stocking of trees > 5 cm DBH, no single definitive standard exists for stems < 5 cm DBH necessary to maintain canopy ingrowth under single-tree selection of northern hardwoods. Older studies frequently evaluated stocking of regeneration irrespective of height class (Leak & Wilson Jr., 1958; Metzger & Tubbs, 1971) or as a total density of stems weighted by height class (McWilliams et al., 1995). 36 Recognizing deer browsing as a formidable barrier to recruitment in some areas, size class must be included when analyzing regeneration success (Walters et al. 2020a). Furthermore, applying standards that were generally developed for even-aged forests, is likely inappropriate for uneven- aged managed northern hardwood stands. Even-aged stocking guides suggest 2850 stems ha-1 for saplings 25 – 137 cm tall, and 950 stems ha-1 for saplings > 137cm tall to < 5 cm DBH (Elenitsky et. al. 2020), or 2500 stems ha-1 for saplings 1.5 - 9.5 cm DBH (Arseneault et al., 2011). Adequate uneven-aged stocking for Arbogast’s smallest size class (~ 4-11 cm DBH) is 494 stems ha-1 after harvest. Within a range of likely stocking guidelines, managers use varying criteria based on harvesting timetables and site-specific factors. When using density of sugar maple saplings to evaluate stocking levels, we found that our stands were largely understocked across conservative criteria (Figure 2.10). For example, based on criteria in Elenitsky et al. (2020), only 29% of stands had > 2500 sugar maple small sapling stems ha-1 and 35% of stands had >1000 large sapling stems ha-1 (Figure 2.10). For a majority (> 50%) of stands to be considered adequately stocked in sugar maple regeneration, thresholds would need to be ~700 stems ha-1 for small saplings and ~500 stems ha-1 for large saplings, the latter is equivalent to Arbogast’s guideline for density of total stems 4 -11 cm DBH. In a silviculture system intended to produce ample sugar maple regeneration, single-tree selection silviculture appears to be failing for many areas in Michigan. It is important to note that since our study did not incorporate unmanaged forests as controls, our results do not necessarily conclude that single-tree selection stands are faring worse than unmanaged northern hardwood stands or that single-tree selection alone has driven these issues. Our finding of widespread understocking differs from Vickers et al. (2019) which found that northern hardwood forests in the Great Lakes region were well stocked. Vickers et al. (2019) included all tree species to evaluate stocking 37 thresholds, including American beech (which cannot presently recruit in high numbers to the canopy due to beech bark disease). Figure 2.10. Percentage of surveyed stands which would be considered stocked in small (50-137 cm tall) and large (137 cm tall to 5 cm DBH) sugar maple saplings, based on average stand densities, is plotted as a function of stocking criteria thresholds (100% stocking at a threshold of 0 stems ha-1, not shown). Red dashed lines provide an example of figure interpretation: at stocking criteria of 2500 stems ha-1 for small sugar maple saplings, 29% of surveyed stands are stocked, and at 1000 stems ha-1 for large sugar maple saplings, 35% of stands are stocked. Given pervasive understocking of sugar maple sapling regeneration and disparate densities of seedlings vs. saplings, it is unsurprising that relative importance of sugar maple saplings is lower than in seedling (Table 2.2) or canopy tree (Table 2.1) classes. Our count data furthermore do not distinguish stems based on vigor or age; old, suppressed regeneration may never recruit to the canopy, even with release from competition (Donoso et al., 2000), artificially inflating sugar maple regeneration density or importance. This may be most pronounced in our large sapling class where many of the stems appeared old and not vigorous (C. Henry personal observation). Altered patterns of relative abundance can lead to compositional shifts (Vickers et al., 2019) which has significant impacts on ecological structure and economic value of northern hardwood forests. Sugar maple replacement with non-viable canopy species, such as beech bark disease infected American beech saplings and emerald ash borer impacted white ash, or by sub-canopy species, such as ironwood, could portend changes in forest structure and loss of a closed canopy 38 (Bohn & Nyland, 2003; Matonis et al., 2011; Bannon et al., 2015; Danyagri et al., 2019; Elenitsky et al., 2020). American beech saplings have particularly been demonstrated to impact understory structure in our region (Elenitsky et al., 2020) and, as one of the most common species in our surveyed stands, may pose future challenges to management. These trends represent significant management challenges in many areas should single-tree selection be maintained as the de facto management system in northern hardwood forests and these trends continue. Our model results highlighted the importance of deer use, site quality, light limitation, and seed availability on sugar maple regeneration and recruitment, and there are several potential management solutions. First, more intensive canopy harvesting regimes that increase understory light availability may improve sugar maple (as well as overall tree diversity) regeneration outcomes (Kern et al., 2017; Webster et al., 2018). Silvicultural systems to potentially achieve this outcome include even-aged shelterwood and seed tree systems, or uneven-aged systems that emphasize larger group selection openings and patch-cuts (Sage et al., 2003; Walters et al., 2016; Hupperts et al., 2020; Walters et al., 2020a). Notably, however, more open overstories are likely to promote greater shrub densities (Walters et al., 2016) and growth of undesirable (for management) advance regeneration (e.g., beech, ash) which compete with sugar maple seedlings and saplings. Thus, it may be necessary to combine more intense overstory harvests with understory treatments aimed at controlling shrubs and undesirable regeneration. In addition to novel silviculture in these systems, treatments may need to include tactics that decrease deer use or browsing pressure on regenerating trees. Although direct management of the deer herd at large scales is difficult to implement and often socially unacceptable, physical barriers, such as felled treetops, may protect regeneration from deer browse (Grisez, 1960; 39 Pellerin et al., 2010; Hagge et al., 2019). Our research also suggests that maintaining larger, seed-bearing sugar maple trees in the stand also positively influences stocking of sugar maple regeneration. Individual managers may select different criterion for designating successful regeneration, but a wide range of thresholds for sugar maple sapling stocking indicates widespread regeneration failure in Michigan. We recommend managers critically examine regeneration in multiple height classes and consider underlying regional drivers to assess potential long-term trends of forest composition and structure; for many regions and forests, continuing with business as usual may lead to fundamental ecosystem shifts. 2.7 Acknowledgements This research was funded by the Michigan Department of Natural Resources. We gratefully acknowledge additional financial and/or in-kind support from Michigan State University, Michigan State University Plant Science Fellowship, Hancock Timber Resource Group, and The Rohatyn Group. The Michigan Department of Natural Resources was involved collaboratively in the study design. Our sponsors were not directly involved in the collection, analysis, or interpretation of data; in the writing of the report; or in the decision to submit the article for publication. We thank our summer interns for their efforts in collecting field data. This chapter was originally published in Forest Ecology and Management (September 2020). The original publication is available at https://doi.org/10.1016/j.foreco.2020.118541. 40 CHAPTER 3 SUGAR MAPLE AGE STRUCTURE IN NORTHERN HARDWOOD STANDS MANAGED BY SELECTION SILVICULTURE 3.1 Abstract Single-tree selection (STS) silviculture has dominated management over the last six decades in northern hardwood forests of the Great Lakes region. Over time, periodic partial harvests are assumed to promote well stocked natural regeneration, resulting in balanced uneven age stand structure with sustainable harvest volumes. However, recent studies indicate understocked tree sapling size classes for species desirable for management, including dominant, economically valuable sugar maple (Acer saccharum) suggests STS in not working, at least recently, in some regions. However, few studies have analyzed sugar maple size-age structure across all size classes to determine whether STS has chronically failed to recruit new age cohorts over the last 60 years or is a recent event. Here, we analyzed stand-scale age structure for 51 managed NHF stands located in the Upper and Lower Peninsulas of Michigan via aging basal discs of 1499 sugar maple trees > 5 cm diameter at breast height (DBH). Our results indicate most stand canopies are dominated by 90 - 120-year-old stems, with low densities of relatively old (> 60 years) saplings (stems 5 – 10 cm DBH). Most cohorts of trees < 90 years old are understocked by uneven-aged stocking guides. For most stands, age and diameter are non-linearly related with a less diverse age structure than assumed when using diameter as a proxy for age. Among stands, by traditional stocking diameter classes, Arbogast’s sapling size class (5 – 11.4 cm DBH) averages 66 years old, poletimber (11.4 – 24 cm DBH) 91 years, and sawtimber (> 24 cm DBH) 106 years. Areas with greater January precipitation (snowfall) and lower annual temperatures, where deer populations have been historically low, generally have younger saplings. Our results 41 offer little evidence to suggest STS has promoted well stocked ingrowth of sapling recruits over 60 + years, i.e., current STS management is unsustainable in the long term in some regions. Our results add a new dimension to the growing body of literature highlighting the need for alternative management of northern hardwood forests in the Great Lakes region. 3.2 Introduction Single-tree selection (STS) has been the predominant management regime for northern hardwood forests (NHF) in the Great Lakes region for nearly six decades (Arbogast 1957, Kern et al., 2014). STS harvests remove individual trees from all diameter classes every 10 – 20 years (Neumann, 2015), creating mostly small harvest gaps (e.g., the width of a dominant tree’s canopy) with the regimen expected to maintain or create (from forests initially even-aged) then maintain a balanced uneven-aged structure. Partial harvesting generates low-light understory environments (Beaudet et al., 2004), which are intended to promote economically valuable (Linehan & Jacobson, 2005; Duval et al., 2014) and ecologically dominant (Schulte et al., 2007; Walters et al., 2020b; Henry et al., 2021) sugar maple (Acer saccharum Marsh) and other shade- tolerant tree species (Crow et al., 2002; Angers et al., 2005; Poznanovic et al., 2013). Despite long-term research and development (Eyre & Zillgitt, 1953; Metzger & Tubbs, 1971), current tree sapling-class patterns in NHF suggest STS may be failing to recruit well stocked stems of desirable species, particularly sugar maple (Leak, 2006; Neuendorff et al., 2007; Powers & Nagel, 2009; Matonis et al., 2011; Henry et al., 2021). What is not known is if sapling-class regeneration failure is a recent phenomenon or indicative of longer-term failure of STS to develop new cohorts following partial harvests. A means of probing this question is to compare existing diameter structures with those expected under uneven-aged STS management and assume that age and diameter are related. 42 However, while diameter structure is important (especially as the currency of stocking guides), shade tolerant sugar maple can persist for decades as suppressed, small diameter stems (Gravel et al., 2011), such that age may vary little among size classes. Employing diameter as a proxy for age may therefore lead to overly optimistic assessments of the uneven-aged character of managed stands. Assessing both diameter and age, and their relationship across size classes, is key for assessing validity of size-based structural guides for STS and the efficacy of STS in promoting regeneration cohorts and recruitment among size/age classes. Here, we use regeneration to refer to the process of establishing a new tree cohort in response to partial harvesting, recruitment as the process of trees transitioning to larger size classes and/or life history stages, and cohort to characterize the group of seedlings/saplings which move into stems > 5 cm DBH per harvesting event. Few studies have directly analyzed sugar maple age structure (Table 3.1), and most are limited in geographic scope and/or diversity of stand structures and underlying landscape-scale drivers. Such analyses, conducted over large regions and considering factors potentially driving pattern, are necessary to determine the extent and reasons for (or associations with) STS successes and failures. 43 Stand history Age at 10 cm DBH (cm) Sampling Location Age-Diam. Source DBH nsamples range height (cm) nstands R2 Tubbs Michigan, USA OG 5 – 81 - 1 60 C 100 N/A (1977a) Leak (1985) New OG 4 – 55 S 1 47 Q 50 0.47 Hampshire, USA Goldblum and Ontario, Canada OG 5-* 30 8 * * * * Rigg (2002) Tubbs Michigan, USA SEL 2.54 – 63 30 1 60 L 48 0.94 (1977a) Kenefic & New York, PAR 2 – 74 137 1 96 Q 27 0.81 Nyland USA T (1999) Dey et al. Wisconsin, SEL 35 - 85 S 4 60 L N/A 0.64 (2017) USA Harmala Michigan, USA VAR > 30 cm tall - 137 8 50 L 42 - 0.82 - (2021) 12.7 cm DBH 70 0.99 Harmala Michigan, USA UM > 30 cm tall - 137 1 5 L 79 0.86 (2021) 12.7 cm DBH Odom and New York, * * 137 * 61 L 24 0.42 Ford (2021) USA Table 3.1. Summary of published sugar maple age-diameter analyses, including source, location, stand history, sampled DBH range, sampling height, the number of stands (nstands), the number of trees sampled (nsamples), the form of the age-diameter relationship, estimated age of 10 cm DBH stems, and the reported R2. Stand history codes include old growth (OG), selection management (SEL), partial cuttings (PART), variable histories, including SEL and PART (VAR), and unmanaged (UM). For sampling height, S refers to stump height (unspecified). Age-diameter relationships are linear (L), quadratic (Q), or free-hand curve (C, which is not a statistical model and therefore does not have an R2 (N/A)). Tubbs (1977a) contained two snapshots of the same stand, before and after harvest treatment, and so is included twice. Harmala (2021) reported stand models separately; for brevity, we have summarized the unmanaged control from the harvest-managed stands, but additional details can be found in Harmala (2021). Information not reported is marked with a *, and information not applicable is marked N/A. 44 The age structure of sugar maple under successful STS management in the Great Lakes region is expected to be influenced by stand history. This includes stand history prior to STS establishment plus the development of subsequent sugar maple cohorts in response to multiple partial harvests over 60+ years of active STS management. Many current NHF in the region established, likely as even-aged stands, 100+ years ago following widespread and intense exploitative logging in the 19th and early 20th centuries (Dickmann & Leefers, 2016), sometimes followed by slash fueled fires (Whitney, 1987). Intense disturbances promoted sprouting species, such as maple, and also early successional species (Whitney, 1987), with the latter since declining in relative abundance due to short lifespans, species-targeted partial harvests, insects and disease (bronze birch borer (Agrilus anxius Gory) impacting paper birch (Betula papyrifera Marsh.), emerald ash borer (Agrilus planipennis Fairmaire) impacting white ash (Fraxinus americana L.)), and limited high intensity disturbance following initial stand establishment (Marquis, 1967; White & Mladenoff, 1994). For most state-owned NHF stands managed by the Michigan Department of Natural Resources (MDNR), partial harvests began in the 1960s with subsequent harvests every 15-20 years. Most partial harvest were STS-motivated, plus many stands in the 1970s were subject to timber stand improvement partial harvests, with similar residual basal area to STS harvests (i.e., 70-90 ft2 acre-1) (Bernie Hubbard, MDNR retired, personal communication). Since the 1960s, most stands likely had 3-4 partial harvests to present. Anticipated age structure of harvest-established cohorts is influenced by sugar maple silvics. Periodic production of large seed crops (2 – 6-year intervals; Garrett & Graber, 1995; Cleavitt & Fahey, 2017) determines initial age distribution of seedlings, which can then persist as seedlings (and eventually saplings) for decades (Marks & Gardescu, 1998; Gravel et al., 2011). However, vigor declines with age/length of low light suppression, such that younger stems 45 demonstrate greater growth in response to harvest than older stems of a comparable size (Donoso et al., 2000), and seedling and sapling populations decline in abundance via accumulated low light mortality slowly with age (Hett, 1971; Gravel et al., 2011). In STS stands, sugar maple can grow to 5 cm DBH in ~ 25 – 40 years (Tubbs, 1977a; Harmala, 2021). Taken as a whole, seedlings-saplings responding to partial harvesting which recruit to a new cohort of stems > 5 cm DBH could include a relatively broad range of ages (several decades) but would likely be numerically dominated by younger individuals from the seedling/sapling class (~ 30 years old). Thus, successful STS in NHF, over the past 60+ years, would be characterized by 3 to 4 harvest-established cohorts, perhaps of relatively broad ages, plus a dominant canopy ~100 - 120 years old. Cohorts of trees recruiting into > 5 cm diameter sapling classes in the last 60 years might be expected to range from (according to Tubbs’ age size relationship) 60-year-old, 15 cm DBH trees to 30-year-old, 5 cm DBH trees. Given sugar maple’s silvics (i.e., as described previously), each partial harvest cohort may vary in age but would likely have a detectable density age peak. Evidence of such age cohort peaks have even been identified in unmanaged NHF, suggesting pulses of recruitment are typical in forests under both managed and unmanaged dynamics (Goldblum & Rigg, 2002). Based on the results in Tubbs (1977a), STS should increase growth rates of poletimber and sawtimber, as evidenced by reducing average age of stems across all size classes; for example, age of a 30 cm DBH stem dropped from 194 years to 104 years, following STS management intervention. Contrastingly, failed regeneration recruitment over the same timeframe would result in a dominant age cohort of stems ~100 – 120 years old occupying a wide range of diameters, from suppressed to dominant stems, with limited younger cohorts. Several factors could contribute to persistent, long-term recruitment failure. Browsing by white-tailed deer (Odocoileus virginianus) could contribute to regeneration failure, with deer 46 populations and browsing pressure varying spatially and temporally and regeneration recruitment failure possibly reflecting these patterns. Temporally, deer populations in Michigan were only ~ 45,000 statewide in 1914, increased to a peak of ~ 1.5 million in the late 1940s, declined to 0.5 million by 1972, increased to a new peak of 2.2 million in 1995, and have declined somewhat but remained relatively high ever since (MDNR, 2016). Thus, we might expect that, with the exception of several years around 1972, deer browsing pressure on young/small trees has been high since the 1940s. Studies have confirmed that high deer browsing pressure is associated with low sugar maple sapling densities (Curtis and Rushmore, 1958; Horsley et al., 2003; Kain et al., 2011; Matonis et al. 2011; White, 2012, Henry et al. 2021) suggesting that high long-term deer densities could chronically limit the recruitment of saplings beyond their reach. In addition to State-scale changes in deer populations over time, deer populations also vary spatially at both local and regional scales due to several factors, some of which vary temporally (Shi et al.2006; Beyer et al., 2010). In Michigan, there is a strong and consistent pattern of low and seasonally migrating deer populations in regions characterized by a deep extended snowpack close to Lake Superior (Beyer et al., 2010). Areas with consistently low deer populations might, all else equal, be expected to have regeneration recruitment less negatively impacted by deer. In addition to deer browsing, several other factors may contribute to persistent sugar maple recruitment failure. The low-intensity harvests of STS may limit recruitment of seedlings to larger classes. STS results in only modest and ephemeral increases in light, with residual canopy tree crowns quickly filling in the small, single-tree canopy gaps and stalling seedlings/sapling growth (Caspersen & Saprunoff, 2005; Kern et al., 2013). Recruitment can also be hampered by light competition from shrubs (Royo and Carson, 2006; Walters et al., 2020a), sedge (Randall & Walters, 2019), and non-sugar maple sapling layers from species like O. 47 virginiana and F. grandifolia (Elenitsky et al., 2020). Finally, although not likely a factor causing failure, site quality as driven by nutrient/water availability may affect the density of sugar maple seedling/sapling populations and the density and composition of its competitors (Elenitsky et al., 2020; Henry et al., 2021). Seedling and canopy tree growth rates are also influenced by site quality, with sugar maple generally having higher growth rate on better quality sites (Kobe, 2006; Baral et al., 2016). Comparing a low vs. high quality NHF, we would expect older saplings (e.g., 5 – 10 cm DBH) and fewer identifiable younger cohorts (i.e., since STS commenced 60 years ago) on the lower quality site due to slower growth rates. Assessing associations of these factors with established sugar maple age structure is key in assessing efficacy of STS management across the range of deer browsing intensities, site qualities/forest community compositions that support NHF. Failure to recruit young age cohorts may portend a future decline in sugar maple canopy stem quality and volume. It has already been reported that poletimber and sapling classes are understocked in Michigan NHF (Walters et al., 2020a; Walters et al., 2022); if existing poletimber and sapling classes are nearly as old as the canopy sawtimber cohort due to limited partial harvest recruitment (i.e., 100-year-old poletimber), then productivity could decline if STS harvesting continues given the preponderance of residual trees with a long history of suppression. In the longer term, with continued STS partial harvests and no recruitment from smaller classes, stocking of the overstory will diminish and ultimately managers will run out of canopy trees to harvest. Exacerbating pole, sapling and ultimately saw timber stocking shortfalls and declining productivity, sugar maple wood quality is likely to decline for sawtimber increasingly recruited from aging, suppressed poletimber and sapling classes (Donoso et al., 2000; Baral et al., 2016; Dey et al., 2017). This would have serious implications for wood 48 product productivity as well as ecosystem functioning and stability. Since stands are estimated to be approaching or surpassing 100 years since stand establishment, and therefore economic maturity for the original age cohort (Dey et al., 2017), the question of success of STS in recruiting younger age cohorts is particularly timely. To address this considerable knowledge gap in our understanding of long-term STS efficacy over the past ~ 60 years of management, we analyzed sugar maple age structure for 51 NHF distributed over a broad area varying in site characteristics, deer use history, and other factors. Given knowledge of management history of forests and deer and current structure of NHF, we predict that single tree selection has failed to recruit new, vigorous age cohorts over the past 60+ years of management. Specifically, we examine the following hypotheses: 1) Age is poorly predicted by typical proxies (e.g., size measurements) a. Age and diameter will be weakly or non-linearly related, with a wide range of tree diameters having similar ages. b. Age will be weakly related to crown class (i.e., suppressed to dominant), live crown ratio, and tree height. 2) Stands will have a single dominant canopy cohort ~ 100 - 120 years old coinciding with the intense exploitative harvests that occurred in the early 20th Century. Except for a possible 45 – 50-year-old class corresponding to low deer populations in the 1970s, there will be little evidence of sapling cohorts originating since partial harvesting began in the 1960’s (e.g., < 60 years old) 3) Observed age structure will differ from theoretical expectations, as evidenced by few trees younger than the age of stand initiation (e.g., 100 years). 49 4) Older saplings will associate with areas of historic high deer density and/or lower site quality. We test these hypotheses by examining sugar maple age-size structure for 51 STS managed NHF stands using aged basal discs from 1499 recently harvested sugar maple trees > 5 cm DBH. 3.3 Methods 3.3.1 Sites and study area To characterize sugar maple age structure, we aged basal disc samples collected from 51 managed NHF stands distributed throughout northern Michigan. Stands were generally on upland mesic to wet mesic sites, with fertile medium-textured upland soils (Dickmann & Leefers, 2016) and dominated by sugar maple on most stands (for total basal area stems > 5 cm, sugar maple most abundant on 48 sites; red maple (Acer rubrum) 2; and balsam fir 1). Stands were considered ready for partial harvest by single-tree selection criteria (~ > 23 m2 ha-1 BA and well stocked in sawtimber classes (Arbogast, 1957). These stands were part of a larger project of 141 stands assessing regeneration outcomes to silvicultural treatments, including four harvesting methods (Walters et al., 2020b) implemented Fall 2017 through Winter 2017/18. This study includes stands harvested by two of the four harvest systems that provide concentrated areas of stumps for sampling: large group selection/patch-cut (numerous dispersed 0.1 – 0.4 ha harvest openings) and seed tree (15 – 20 tree ha-1 retained; Walters et al. 2020b). From a pool of 71 seed tree and large group selection/patch-cut harvested stands, we stratified sampling by three geographic areas (Western Upper Peninsula, Eastern Upper Peninsula, Northern Lower Peninsula, n=18 per region) and three ordinal site qualities (high, medium, and low, categorized by a habitat classification system based on understory flora (Burger & Kotar, 2003)) to select the 54 stands for this study. Three stands were dropped due to harvesting or logistical challenges, 50 resulting in a total of 51 stands. Due to dropped stands and variable representation of stands by site qualities across regions, representation of stands by strata were imbalanced (Appendix D, Table D.1). 3.3.2 Field methods Prior to harvest, we established plots in summer 2017 to characterize stand structure and composition; we also randomly selected sugar maple trees for age analysis to measure pre- harvest individual tree characteristics, such as crown class (i.e., relative competitive status). From a permanent grid of 25 survey points, we randomly selected plots with minimal slope (< 10 degrees) and adequate sugar maple tree sample size (> 10 sugar maple trees, > 10 cm DBH, within a 25 m radius plot), randomly resampling if plots failed these criteria. For large group selection/patch-cut harvest stands, random sampling was restricted to the two largest harvest opening sizes (centered on 10 of the 25 survey points). To characterize diameter distribution and composition of canopy stems directly competing with our sampled trees, we recorded DBH, species, and alive/dead status for all trees > 10 cm DBH within a 25-m radius (0.196 ha) plot (referred to as plot-level data). To characterize the broader diameter distribution and composition of the stand, we utilized a larger dataset of stem DBH, species, and alive/dead status (25 plots dispersed throughout the 12.14 ha stand, each 113 m2, totaling 0.283 ha per stand; hereafter referred to as stand-level data). Since we wanted to sample age equally amongst size classes, we used stratified sampling within three DBH sampling bins of equal width, ranging from 10 cm DBH to the plot-specific second-largest sugar maple DBH, to randomly select 21 sugar maple stems. We also randomly selected up to 7 sugar maple stems 5 – 10 cm DBH to ensure saplings were well represented in the data. We maximized spatial dispersion of sampled trees by selecting the closest stem in each 51 of the four bins at 7 spatially dispersed sub-sampling positions within 20 m of plot center. We attempted to restrict tree selection to within 20 m of plot center to minimize competitive influence of trees just beyond the 25 m radius plot which are not quantified in our diameter distribution data. For each selected stem, we recorded DBH, crown class (dominant, co- dominant, intermediate, overtopped, adapted from Nyland, 2002), and live crown ratio (Schomaker et al., 2007), plus height for intermediate and overtopped trees, and labeled each with an aluminum tag near the base. For the three bins > 10 cm DBH, if fewer than 7 trees were located, samples were re-allocated to other size classes to attempt to always tag 21 trees. In summer 2018, following harvest, we cut cross-sectional discs from tagged stumps. Many tags placed in summer 2017, particularly for smaller stems, were missing following harvest, most likely due to logging equipment impacts. Therefore, we additionally sampled up to 12 untagged stumps, as needed, equally allocated among the three size bins > 10 cm DBH as previously described, excluding stumps with obvious missing piths due to decay. Up to 5 untagged saplings 5 – 10 cm DBH were also sampled, often by searching a larger area including areas outside, but close, to plots if necessary. We generally collected entire basal discs, but in some cases were only able to recover a portion (but still including pith and most recent rings) due to partially shattered stumps, stumps cut low to the ground, and large stumps. To capture general variance in sampling height, we recorded an estimated stand-level average height for stumps. We collected a total of 1521 sugar maple basal discs. 3.3.3 Lab methods We dried and sanded basal discs with progressively finer sandpaper (up to 400 grit) to reveal annual rings. We counted rings on two radii per sample, with the first radius typically the longest radius on the sample, which had most easily identifiable (widest) rings and was least 52 likely to have missing rings (Lorimer et al., 1999). If the ages of the two radii differed by more than 5 years, we counted a third radius to help ensure we identified the maximum number of rings present. We excluded samples with more than 2 cm radius of pith missing (n = 22, mostly due to rot), as pith ring widths vary widely and may not be estimated reliably (personal observation). Our final sample pool was 1499 discs, an average of 29 samples per stand (range 19 – 38). We were unable to locate any saplings (5 – 10 cm DBH) in three stands, including one stand in which we were unable to locate any stems < 20 cm. 3.3.4 Data analysis 3.3.4.1 Estimating age and diameter at breast height Sugar maple is more likely to have missing versus false rings (Lorimer et al., 1999), including for defoliation events, which trigger reduced growth rates as opposed to false rings (Bevilacqua et al., 2021). Therefore, we estimated tree age as the maximum number of rings counted from a single radius. For basal discs missing less than 2 cm radius to the pith (n = 63), we estimated the number of missing rings by estimating distance to the missing center and dividing that distance by the average ring width of the inner-most 20 rings (Duncan, 1989). To ensure that plot-level stump height was not influencing our age, we modeled age as a function of stump height; there was no significant effect (p = 0.61), so we took no further action to correct for variation in stump height. DBH is a universally used field measurement by ecologists and foresters, and which we used to characterize pre-harvest structure data; we therefore wanted to use DBH rather than disc sample diameter for modeling and age-diameter analyses. Since about half (739 of the 1499) of our basal discs were from untagged trees and lacked pre-harvest DBH measurements, we assessed whether basal disc diameter accurately estimated DBH. We measured up to 6 radii 53 inside the bark, evenly spaced at approximately 60° intervals, for all basal discs; fewer radii were measured for incomplete or damaged basal discs. For basal discs from tagged trees, we modeled DBH as a function of average basal disc diameter using the lm function in R (n = 760). Since average basal disc diameter was a strong estimator of DBH (R2= 0.85, yDBH = 1.35 + 0.86 * xdiam, where yDBH is sample DBH in centimeters and xdiam is basal disc diameter in centimeters), we calculated estimated DBH based on basal disc diameter for all samples in place of incomplete field DBH measurements, and heretofore refer to this as sample diameter or DBH interchangeably. 3.3.4.2 Age associations with tree size characteristics (Hypothesis 1) To quantify stand-level age-diameter relationships, we compared a linear regression model (fit with lm function in R) and a non-linear Chapman Richards growth function model (Richards, 1959; D. G. Chapman, 1961). To fit the Chapman Richards model, we used nlsLM from the minpack.lm package (Elzhov et al., 2016), which incorporates the Levenberg- Marquardt fitting algorithm (Levenberg, 1944; Marquardt, 1963), and generated model starting values with nls_table from the forestmangr package (Braga et al., 2021). We used Akaike's Information Criterion (AIC) (Sakamoto et al., 1986) to compare models. We additionally considered quadratic, exponential, and segmented linear regression, but these models provided biologically implausible fits (quadratic or exponential), were highly influenced by outlier or leverage points (segmented linear regression), or fit poorly (exponential). We visually assessed preliminary models for outliers and high leverage values and removed five samples for final model fitting. To characterize variation explained by chosen models for each stand, we calculated Efron’s pseudo R2 values using the accuracy function from the rcompanion package (Mangiafico, 2021). In addition to individual stand modeling, we pooled samples from all stands 54 to develop an age-diameter relationship for the entire study region, using the same methodology as for individual stands. For all other tree-specific factors assumed predictive of age, we ran three separate mixed effects models of age using a normal distribution and including stand as a random effect. We used lmer from the lme4 package (Bates et al., 2015). Predictors of the three models were: canopy (4 levels); live crown ratio (expressed as a percentage) and DBH, with an interaction effect to account for size-related trends in live crown ratio; and height. The first two models, respectively, were conducted for the full set of tagged and recovered samples (n = 760), while the third model included only overtopped and intermediate tagged and recovered stems with available height data (n = 300). We calculated marginal R2 values (e.g., the variance explained only by the fixed effects) for all models using r.squaredGLMM from the MuMIn package (Barton, 2020). To summarize age structure in a format directly applicable to management, we summarized median stand age of sampled discs in each of Arbogast’s timber size classes (subset of Arbogast saplings which we sampled, 5 – 11.4 cm DBH, plus poletimber (11.4 – 24.13 cm DBH) and sawtimber (> 24.13 cm DBH)) for each stand. 3.3.4.3 Age cohort analysis (Hypothesis 2) To characterize age cohort density, we generated expected stand-level density distributions (stems ha-1) for sugar maple age classes. We used stand-level, rather than plot-level, data to characterize expected age density distribution due to greater area surveyed in stand-level data and to make the results applicable at a stand-scale. We assigned each basal disc a representative stem density based on stand-level diameter distribution data, using 5 cm diameter bin widths to crosswalk basal disc samples with sugar maple stand diameter distribution. For 55 example, in a given stand, if the stand-level density of sugar maple trees 10 – 15 cm DBH was 20 stems ha-1, and we sampled 2 basal discs 10 – 15 cm DBH, each basal disc is assigned a representative density of 10 stems ha-1. Given sparsity and skew of large diameter trees, we defined diameter groups in 5 cm intervals from 5 cm DBH to the diameter of the second largest sugar maple basal disc sample, rounded down to the nearest 5, and made the final diameter bin open-ended to include all large stems (for both basal discs and plot-level diameter structure). We occasionally combined adjacent diameter bins when we lacked basal discs in a given diameter class. Expected age densities were then aggregated by age bins. This method may slightly underestimate variability of ages present in the stand since it relies entirely on the age of collected discs, particularly in cases where diameter classes were under-sampled (e.g., smaller stems more likely to be fractured by logging equipment). To identify potential age cohorts, we characterized peaks in the age-class density distribution of our samples. We generated five histogram iterations by shifting bin centers by 1- year increments using 5-year age bin widths, since histograms are sensitive to both bin width and position (Pond & Froese, 2015). Repeated for each histogram iteration, an age-class bin had to 𝐷 surpass a threshold density to count as a peak. Our threshold density, T, was defined as 𝑇 = , 𝑛 where D is the total sugar maple stem density and n is the number of age bins (from minimum to maximum age bin, including unoccupied age bins). We developed T via an iterative process, testing additional thresholds of 10% of total sugar maple basal area and average bin density excluding empty bins. T successfully characterized peaks in our data without being overly sensitive, and peaks were generally robust to varying threshold criteria. To determine the width of an age cohort, we considered adjacent bins greater than T, plus 1 any adjacent bins greater than one half the value of the threshold ( 2T), as a single age cohort. 56 1 Including lower density bins (2T) to define cohort width accounts for the anticipated diffuse nature of age cohorts due to competition during stand establishment and due to variability in aging (e.g., missing rings), sample size, and stump height. Histogram iterations were generally robust in identifying similar peaks and age cohorts. Therefore, for each individual age cohort, we defined its placement and width by selecting the histogram iteration which captured the peak in the least number of age bins, randomly selecting when iterations tied for narrowest width; in doing so, we maximized the number of age cohorts which could be defined with our criteria. After age cohorts were defined for each stand, we visually identified one age cohort which captured most co/dominant stems (typically > 20 cm DBH) which we hereafter refer to as dominant canopy age. We calculated Moran’s I to assess landscape spatial autocorrelation in dominant canopy age among stands (Moran, 1950). 3.3.4.4 Comparison with theoretical age structure (Hypothesis 3) To develop expected stocking by age class (as described in 2.4.4), we combined density stocking guidelines outlined in Arbogast (1957), a widely employed stocking guidance manual, with the sugar maple age-diameter relationship quantified in Tubbs (1977a), one of the seminal works establishing age-diameter assumptions for sugar maple in STS managed NHF. Using the equation in Tubbs (1977a), we estimated age for stems ranging from 5 – 61 cm DBH (DBH range of Arbogast guide which we sampled) at 0.25 cm intervals, equivalent to sampling ten stems within each Arbogast 1-inch size class. We assigned each stem an expected density using Arbogast’s estimated stem density per hectare for 1-inch size classes, assigning each stem as 1/10 of the total stem density of that class; note, we did not fully sample the smallest Arbogast size class (1.5 – 2.5 inches DBH), thus we are excluding stems < 5 cm from our Arbogast theoretical expectations. We then aggregated stems ha-1 for age bins to match our data 57 categorizations. Although Arbogast’s guidelines are intended for all stems and not exclusively sugar maple, past management emphasis for sugar maple and dominance of sugar maple in current stands, plus the fact that Arbogast’s guidelines are considered minimum residual stem densities targets following harvest (and we are using pre-harvest data), makes this a reasonable approximation of stocking densities and an accurate description of expected age distribution shape. 3.3.4.5 Stand and landscape factors associated stand age structure (Hypothesis 4) We wanted to explore how stem age for key diameter classes varied by stand and landscape-level predictors. To maximize comparability with other studies and minimize sampling bias, we extracted model predicted age values for stems which are 10 cm, 20 cm, and 30 cm DBH (referred to as AGE10, AGE20, and AGE30), per stand. These three DBHs represent key diameter classes: large, recruited saplings (AGE10), poletimber (AGE20), and sawtimber (AGE 30), and were generally within the range of diameters sampled per stand. We modeled AGE10, AGE20, and AGE30 as a function of two long-term climate variables (PRISM 30-year normal estimates (1981 – 2010) of annual mean temperature and January precipitation (snow); PRISM Climate Group, Oregon State University, http://prism.oregonstate.edu, created 23 September 2021) and two current stand metrics: plot- level tree density (total basal area of stems > 10cm) and site quality (ordinal variable, ranked 1 – 4 with 1 lowest and 4 highest, Burger and Kotar 2004; Table 3.2). Since the Great Lakes generate long-term (e.g., centuries) predictable regions of deeper snow (i.e., lake effect snow) than surrounding areas and because winter temperature and snow are well-established deer population drivers (Moen, 1976; Holter et al., 1975; Shi et al., 2006), we use normalized climatic measures as a proxy for historic deer population distribution. Current stand-level characteristics 58 reflect the culmination of management and other disturbances/effects and growth over the last several decades. Site quality indicates relative differences in growing conditions over much longer time scales, driven largely by soil textural differences among post-glacial landforms (Zak et al., 1989; Baribault & Kobe, 2011). Variable Mean (Range) / n Annual mean temperature 5.6 (4.6 – 7.6) °C January precipitation 51.8 (29.4 – 74.5) mm Total basal area of stems > 10 cm 26.8 (14.8 – 54.0) m2 ha-1 Site quality 1 n = 15 Site quality 2 n=9 Site quality 3 n = 17 Site quality 4 n = 10 Table 3.2. For predictors of AGE10, AGE20, and AGE30, mean and range of for continuous variables and number for categorical site quality. Site quality 1 is lowest and site quality 4 is highest. We calculated Moran’s I (Moran, 1950) to assess for spatial autocorrelation; only AGE10 was spatially autocorrelated (observed = 0.092, p-value = 0.029). We implemented regression modeling using a normal distribution with the spLM function from the spBayes package, which uses a Bayesian framework and can incorporate a spatial term if needed (Finley et al., 2015; Finley & Banerjee, 2020). We standardized our continuous predictor variables and ran models for 25,000 iterations. To check model convergence, we assessed chains for evidence of mixing. We calculated Moran’s I on residuals for the AGE10 model; the model was able to account for all spatial autocorrelation, so we did not include a spatial term in our AGE10 model. 59 3.4 Results 3.4.1 Data overview Local growing conditions within plots varied in basal area, sugar maple dominance, and percentage of standing dead trees (Table 3.2). Sugar maple density ranged from 102 to 652 stems ha-1 with relative abundance ranging from 16 – 100% (average 72%; Table 3.3). We tallied 20 total tree species across the 51 plots/stands (Appendix D, Table D.2). Based on stand-level data, diameter distribution was generally understocked (Arbogast 1957) in stems 5 – 10 cm DBH and stocked to overstocked in stems 10 - 45 cm DBH (Figure 3.1). Figure 3.1. Stand-level diameter distribution of sugar maple, conifer, undesirable, and other desirable species in 5 cm diameter size classes (size classes are labeled with median diameter value) from 51 northern hardwood forests stands sampled in northern Michigan, USA. Undesirable includes American beech and ironwood, while other desirable refers to all species other than sugar maple and undesirable species (sensu Walters et al., 2022 but excluding ash as undesirable species). Black asterisks represent theoretical diameter class stem densities at full stocking, for all species > 5 cm DBH, according to Arbogast (1957). 60 Variable Mean (range) Total stem density 420 (163 – 917) stems ha-1 Total BA 27 (15 – 54) m2 ha-1 Sugar maple stem density 295 (102 – 652) stems ha-1 Sugar maple BA 18 (3 – 54) m2 ha-1 Sugar maple relative abundance, by stem count 72 (16 – 100) % Sugar maple relative abundance, by BA 68 (14 – 100) % Table 3.3. Plot-level diameter composition of all living stems > 10 cm DBH for 51 plots, each representative of stand, in which age discs were sampled. Ages ranged from 23 to 255 years, with a median and mean of 96 years. The youngest stem sampled in each stand averaged 52 years (range 23 to 118 years). The smallest stem sampled in each stand averaged 7.4 cm DBH (median 6.5, range from 5.5 – 28.4; saplings were present in our stand structure data but not located for age sampling in three stands, including one stand in which no stems < 20 cm DBH were sampled). 3.4.2 Associations between age and size characteristics (Hypothesis 1) For data pooled across all stands, the Chapman Richards function was selected over a linear model (AIC 13461 vs. 13503, respectively) for fitting age vs. diameter, with slope decreasing as diameter increased (e.g., a wide range of larger diameter stems with similar ages) (Figure 3.2). For individual stands, AIC scores supported Chapman Richards model fits for 36 stands and linear model fits for 15 stands (Figure 3.3; Appendix D, Table D.3). For the 36 stands modeled by the Chapman-Richards equation, age plateaued at 111 years, on average (range from 93 - 157 years) (Appendix D, Table D.4). For the 15 stands modeled by a linear equation, two slopes were not significantly different from zero (p > 0.05), and significant slopes averaged 2.1 (range from 0.4 to 4.3); all but one intercept were significantly different from zero, and significant intercepts averaged 58 years (range from 31 – 118 years) (Appendix D, Table D.5). 61 Pseudo r-squared values varied widely for all models, ranging from 0.00 to 0.92 with a median of 0.63 and mean 0.58 (Appendix D, Table D.3). Figure 3.2. Age versus diameter (cm) for 1499 sugar maple stems ≥ 5 cm in managed northern hardwood forests (black points). The Chapman Richards model fit to our data is shown in orange. The other lines represent relationships between age and diameter for sugar maple stems from other studies (for the range of diameter values in each study). Details on studies can be found in Table 3.1. 62 Figure 3.3. Model-predicted values of sugar maple age by diameter (for the stand-specific range of diameters) for 51 managed northern hardwood forests. Black lines are linear models and red lines are Chapman-Richards models. For all data (i.e., stands pooled), age varied among crown classes, with means of 122, 105, 94, and 80 years for stems in dominant, codominant, intermediate, and overtopped (i.e., suppressed) positions respectively (marginal R2 = 0.22; Appendix D, Table D.6). There was considerable variation in age within each crown class except the overtopped class (Appendix D, Figure D.1). Tree diameter and live crown ratio interacted to predict stem age. For stems with smaller diameters, age tended to decrease with increasing crown percentage, and for stems with larger diameters, age tended to increase with increasing live crown ratio; however, the model (age as a function of diameter, live crown ratio, and their interaction) was relatively weak (marginal R2 = 0.32; Appendix D, Table D.7). Age also positively associated with height (estimated model, A ~ 55.7 + 2.2*Hm, where A is age, Hm is height (m)), although weakly (margin R2 = 0.22; Appendix D, Table D.8). Among stands, median age of sampled stems in Arbogast’s three size classes (saplings, poletimber, and sawtimber) varied (Figure 3.4). The average median value of sapling discs was 69 years, with most stand median values between 50 – 100 years but eight stands (of 49 with sapling age samples) had median sapling ages < 50 years. Stand-level median age of poletimber 63 appears bimodal, with peaks around ~85 and ~100 years (overall mean 90 years). Sawtimber had a median age of 102 years (mean 106) and ranged from 72 to 166 years. Figure 3.4. Histogram (5-year bins) and boxplot of median age for sampled discs in Arbogast’s sapling (5 – 11.4 cm DBH), poletimber (11.4 – 24.1 cm DBH), and sawtimber classes (> 24.1 cm DBH), by stand, for 51 STS managed northern hardwood forests. Dashed lines emphasize 50, 100, and 150 years on the x axis. Boxplots show first, second (median), and third quartiles, plus whiskers which extend up to approximately three standard deviations from the mean (1.5 times the inter-quartile range); points represent outliers. Note, we did not sample the full width of Arbogast’s sapling class (stems 3.8 – 11.4 cm DBH). 3.4.3 Age cohort results (Hypothesis 2) Age distributions by stand were variable but showed some general trends (Figure 3.5). Most (59%) stands had three or fewer age cohorts (8 stands had one cohort; 14 had two; 8 had three). The median age of the dominant canopy cohort averaged 101 years, ranging from 53 to 163 (Figure 3.5C); the width of the dominant canopy cohort (e.g., number of age bins it spanned) ranged from 5 to 35 years (median 15 years, or three 5-year age bins). There was no statistically 64 detectable spatial autocorrelation in the dominant canopy age (Moran’s I, p value = 0.15, Figure 3.5C). For 12 of 14 stands with two cohorts, there was a cohort ~ 20 – 30 years younger than the main canopy cohort, with all but one of these younger cohorts > 50 years old. Cohorts younger than the canopy were typically low-density, in opposition of J-shaped size-age distribution expectations which with greater stem density in small/younger cohorts. Among all stands, few cohorts (n = 17 of 170, or 10% of all identified cohorts) had a median age < 50 years old. Sixteen stands had at least one cohort older than the dominant canopy class, including five stands with age cohorts above 200 years located in the Upper Peninsula (4 of 5 along the northern half). 65 A B C Figure 3.5. Histogram-identified sugar maple age cohorts by site (x-axis) for 51 managed northern hardwood stands (3.5A). Point size represents total stems per hectare measured within that cohort, and black points represent the age cohort capturing the dominant canopy. Stands are organized by number of identified cohorts within each stand and sorted (low to high) by median age of the youngest age cohort. Number of identified age cohorts (3.5B) and median age of the canopy cohort (3.5C) are represented spatially, plus a histogram of canopy ages (3.5C). 3.4.4 Age structure deviation from expectations (Hypothesis 3) For all stands pooled, stand-level sugar maple stem density by age class appears normally distributed, with a modal age class of 90 – 100 years old (columns, Figure 3.6A). Compared to theoretical age distribution for successful STS management (black dots, Figure 3.6A), age classes < 60 were understocked, age classes and age classes 80 - 120 years old were stocked to overstocked, and age classes > 120 years were understocked (Figure 3.6A). Stand-level analyses 66 indicated greater variability in stocking by age class (Figure 3.6B). The 70 – 90 and 90 – 110- year age bins were most often stocked, but 25 and 22% of stands, respectively, were not stocked in these age classes, though they typically had a surplus of stems in an adjacent age class (Figure 3.6B). No stands were stocked in the 30 – 50-year age bin (Figure 3.6B). A B Figure 3.6. Bar plot of average age density distribution (20-year age bins, labeled with median value) of sugar maple stems > 5 cm DBH, pooled across 51 managed NHF stands; black points represent age density of a fully stocked northern hardwood stand, following harvest, if entirely stocked with sugar maple (Figure 3.6A). Anticipated age distribution was determined by applying the age-diameter relationship from Tubbs (1977a) to the size class densities outlined by Arbogast (1957). Bar plot representing the percentage of stands individually considered stocked (e.g., exceed Arbogast’s extrapolated stocking value) for each age class expected under the anticipated age distribution (Figure 3.6B). 3.4.5 Factors associated with stand age structure (Hypothesis 4) Based on stand level age-diameter models, predicted ages of stems at 10, 20, and 30 cm DBH (AGE10, AGE20, and AGE30) averaged 72, 94, and 106 years, respectively, with ranges of ~50+ years (Figure 3.7). AGE10 stem age was greater in areas with higher mean temperatures and lower January precipitation, i.e., regions with milder winters, where we expect long-term deer browsing pressure to be higher (Table 3.4). AGE20 stem age also increased with lower January precipitation, and while AGE30 was unrelated to climate variables, stems on medium quality sites were predicted to be older than stems on low quality sites (Table 3.3). 67 Figure 3.7. Kernel density estimates of model-estimated values of stems 10 cm (purple solid line), 20 cm (green dashed line), and 30 cm (yellow dotted line) DBH are presented (n = 48 stands for AGE10, n = 50 for AGE20, and n = 51 for AGE30); these are interpretable as a smoothed version of a histogram, where the area under the curve of each line sums to one. We did not extract model estimates beyond the range of data, therefore estimates were not available for three stands lacking saplings, one of which additionally had no stems below 20 cm DBH. 68 Size Variable E 50% E 2.5% E 97.5% Intercept 71.5 64.9 77.7 Mean annual temperature 6.4 2.6 10.4 January precipitation -5.5 -9.4 -1.4 Plot basal area 0.7 -3.1 4.6 AGE10 Site quality 2 0.8 -9.3 11 Site quality 3 1.6 -7.2 11.2 Site quality 4 0.1 -10 10.4 T2 145 97.3 232 Intercept 93.5 87.4 99.5 Mean annual temperature 3.1 -0.5 6.9 January precipitation -3.9 -7.7 -0.3 Plot basal area 0.9 -2.8 4.3 AGE20 Site quality 2 5.3 -4.6 15.1 Site quality 3 -0.3 -9.1 8.4 Site quality 4 -1.8 -11.3 8.3 T2 131 87.7 207 Intercept 102.7 95.5 109.7 Mean annual temperature 2.2 -1.8 6.4 January precipitation -0.1 -4.1 4.1 Plot basal area -1.2 -5.1 2.6 AGE30 Site quality 2 14.9 3.6 26.3 Site quality 3 2.8 -7.2 13.2 Site quality 4 -2.4 -14 9 T2 169.2 114.8 264.4 Table 3.4. Model parameter estimates (E 50%) and 95% credible intervals (E 2.5%, E 97.5%) for factors predicting age of sugar maple stems at 10, 20, and 30 cm DBH. Significant predictors (p < 0.05) are bolded and italicized. Site quality 1 (lowest site quality) is included in the intercept, and site quality 4 is the highest site quality. T2 represents the variance (equal to the squared standard deviation). 69 3.5 Discussion 3.5.1 Overview Age-diameter relationships from 51 STS managed NHF in northern Michigan support our overall hypothesis that STS has largely failed to recruit new regeneration cohorts (> 5 cm DBH) over the past ~ 60 years of management. Successful uneven-aged management is expected to generate linear tree diameter vs. age relationships (such as in Tubbs, 1977a), but our results suggest stems of similar age can occupy a wide range of pole and sawtimber diameter classes. Sapling classes are generally understocked by STS management standards (Arbogast 1957) and old (5 – 10 cm DBH stems average 66 years old). Overall, few sugar maple trees recruited (grown into > 5 cm DBH classes) since STS became widespread, indicating extensive regeneration and recruitment failure under STS management for sugar maple. Our results support Hypothesis 1, that age is poorly estimated by diameter and growth form characteristics (crown class, height of understory trees, and live crown ratio); small diameter stems with small crowns in subordinate crown classes were as, or nearly as old, as dominant, large diameter, and large crowned canopy stems. Diameter – fundamentally the most common age proxy – often displayed a non-linear relationship with age, and for the 15 stands best fit by a linear model, slopes were often very small (and two were not significantly different from zero; Appendix D, Table D.2). Variation among relationships indicates that age-diameter relationships need to be quantified at sub-regional (if not stand-level) scales if they are to be used as an accurate, useful tool for estimating age from diameter. Our second hypothesis that stands are dominated by a canopy cohort of sugar maple stems, approximately 100 – 120 years old, was supported by the data (Figure 3.5A). Interestingly, our results suggest sapling recruitment of current dominant canopy trees (under 70 mostly open, highly disturbed conditions) was typically gradual (e.g., average width 15 years). Extended canopy establishment may be explained by competition from both non-tree vegetation (e.g., Rubus shrubs, Hughes & Fahey, 1991; Walters et al., 2016) and faster growing early successional trees (aspen, birch; Marquis, 1967) delaying/extending establishment and growth of sugar maple. A mixture of stump sprouting and seed-origin establishment, which vary in vertical growth rate and establishment time (Whitney, 1987), could also contribute to this phenomenon. Canopy age was not spatially autocorrelated across our study region, which suggests local factors drive canopy establishment rather than regionally varying factors. Despite lack of statistical evidence, visual inspection of mapped canopy age indicates most stands with older canopy ages (~ 125 years old) are located along the Great Lakes and near areas with long settlement history (e.g., Traverse City, Houghton), while younger canopies tend to be found in more remote locations (Figure 3.5C). This may suggest that older forests are associated with areas with greater accessibility to mills early in the exploitation era and/or areas where earlier conversion to agriculture followed by abandonment occurred (Dickmann & Leefers, 2016). The characteristics of additional age cohorts for stands with more than one age cohort further support Hypothesis 2; there was little evidence of young age cohorts establishing since partial harvesting began under STS in the 1960s. Although some stands had 4+ age cohorts, most were older (often older than main canopy cohort). The common pattern of a cohort 20 - 30 years younger than the main canopy cohort defies easy explanation. This secondary cohort may have filled the growing space initially filled by a short-lived early successional species that subsequently died (e.g., pin cherry; Whitney, 1987). They may also reflect the first harvest of a stand that generated a new age cohort; however, they mostly precede, in time, anecdotal evidence of partial harvests (1960s), nor is there any conceivable justification for partially 71 harvesting a 30-year-old stand. Additional stand history investigation, possibly combining detailed records, if they existed, with age distribution data, may yield greater insight. We were further surprised to identify several stands with sugar maple age cohorts > 200 years old, generally in the northern Upper Peninsula. That those stands were not entirely clearcut during the era of exploitive intense harvests (late 1800s – 1920s) is not surprising given that exploitative harvests of that era may have often resulted in non-commercial (small, poor form) stems being left behind. Alternatively, some stands likely had old-growth structure with no harvesting up to the time that STS was implemented in the 1960s; three of our stands which had 5+ age cohorts of similar size including age cohorts > 200 years are potential candidates for this category (Figure 3.5A). The widespread lack of sugar maple cohorts < 50 years old coincides temporally with increasing deer densities beginning in the 1970’s, related to mill expansion and the creation of abundant browse via widespread harvesting of aspen (MDNR 2016). However, the lack of cohorts 50 – 60 years old defies expectations, as establishment of these stems in the 1960’s would have coincided with relatively low deer populations (MDNR, 2016) and the beginning of STS partial harvesting (Bernie Hubbard, MDNR retired, personal communication). The lack of saplings (e.g., > 5 cm DBH) aged 50 – 60 years suggests STS widely failed to create well stocked sapling recruit classes throughout northern Michigan despite lower deer populations. While deer browsing likely still contributed to low sapling recruitment other factors may have also contributed, including small STS gaps creating insufficient light environments for sugar maple growth and survival via maintenance of high basal area/low light (Henry et al., 2021). In addition, STS commenced when stands were young (40 - 50 years) with likely small diameters, such that sugar maple seed (the source of future sapling recruits) could have been limiting in 72 many stands (Henry et al., 2021). Overall, our results point to a surplus of aging sugar maple canopies and a lack of younger recruited age cohorts, indicative of persistent recruitment failure. Lastly, we found support for our hypothesis that high deer use would be associated with older saplings but not for high site quality associating with younger saplings (Hypothesis 4). AGE10 was spatially autocorrelated, with regions of younger sapling ages along the northern half of the Upper Peninsula where snow is deep in the winter (Appendix D, Figure D.2). Related, AGE10 was older in regions with warmer winters and less snow, with both factors known to be associated, in parts of the study region, with higher winter deer populations in the region (Holter et al., 1975; Moen, 1976; Shi et al., 2006). Given the general lack of AGE vs. site quality relationships, we found no support for our hypothesis that presumably higher growth rates on higher quality sites would lead to younger stems at given diameter (e.g., lower AGE10). This could be due to lack of a site quality vs. growth relationship or other factors (deer browsing, stand density) may have overwhelmed or confounded site effects. As an example of potentially confounding effects, basal area and shading of subordinate tree strata could be greater on higher quality sites and diminish growth. Regardless of factors driving pattern, overall, the preponderance of old, low density sapling classes in managed NHF is concerning and emphasizes the importance of understanding age as well as diameter distributions for stocking if employing STS. 3.5.2 Comparison to other studies Our size-age relationships differ from previous studies of sugar maple in managed forests (Table 3.1) and appear more similar to those for old-growth/unmanaged forests (Tubbs, 1977a; Canham, 1985). Compared to previous studies conducted in the Great Lakes region, differences in our results might be explained by our study encompassing a broader region (as compared to 73 Harmala 2021) with more variation in deer use (as compared to studies in areas of low deer populations Dey et al. 2017 and Harmala 2021). Restricting our sampling to one or two stands, particularly in a region of lower deer browsing, may have resulted in similar conclusions; for example, the map of AGE10, which was spatially autocorrelated, indicates a cluster of stands along the central northern portion of the Upper Peninsula which have stems 10 cm DBH averaging 50 – 60 years old (Appendix D, Figure D.2), which is more comparable to past studies (Table 3.1). Our summary also indicates macro-regional variation in sugar maple age-diameter relationships, with stems older at given size in Michigan (Tubbs, 1977a) than in New England, both for old-growth (e.g., Leak 1985 in New York) and managed stands (e.g., Kenefic and Nyland 1999 in New York). This may be explained by greater summer precipitation in the Northeast (~76 – 89 cm annual precipitation in northern Michigan vs. 102 – 152 cm in New England, PRISM 30-year normals) or by regional variation in stand composition and history, including management. However, contrary to patterns for larger trees, Kobe (1996) reports similar growth rates for sugar maple saplings in Michigan vs. Connecticut, USA. An important caveat to our comparison of studies is that sampling height can influence age estimations, since it can take sugar maple a median of 2 – 10 years to progress from 30 cm to breast height (Harmala, 2021). Odom and Ford (2021) and Kenefic and Nyland (1999) sampled at breast height; however, adding ten years to their estimates of sugar maple age (across all diameters) would still indicate stems in New England are considerably younger than in Michigan. Overall, except for saplings, evidence suggests sugar maple grows more quickly in New England than the Great Lakes region. 74 3.5.3 Management implications Individually and in combination, stand diameter and age distributions are concerning for management by STS (Figures 3.1, 3.6). The deficit of young stems in sapling and poletimber diameter classes suggests persistent recruitment failure, which has several ecological implications. Older stems tend to acquire greater damage and release with less vigor (e.g., slower basal area increment) in response to harvest (Baral et al., 2016), such that stands dominated by old stems could have lower growth rates than stands dominated by young stems. Lack of vigorous young stems may also mean that sugar maple is less likely to capture gaps following canopy disturbance or harvest. If this provides opportunity for species diverse stems to recruit, this outcome could be beneficial for stand resilience; however, it could also lead to gap capture by highly abundant, non-viable (American beech due to beech bark disease) or non-commercial (ironwood) species (Elenitsky et al., 2020; Walters et al., 2022), which would lead to canopy failure (e.g., loss of continuous canopy). From a timber supply standpoint, our results point to pressing challenges in both near (e.g., one to two decades) and intermediate (e.g., four to five decades) time frames. Poletimber stems averaging ~ ten years younger than sawtimber indicates an imminent impact of unbalanced age structure. Given that sugar maple stems tend to decline in quality after 100 years due to decreasing proportions of high-grade white wood lumber (Dey et al., 2017), trees harvested as larger sawlogs (e.g., 40 cm DBH +) in the future will have lower rates of return than those harvested at 40 cm DBH now (i.e., at approximately 100 years old). Stands with limited age distribution may benefit from more intense harvests to release poletimber to the sawtimber class and to increase light availability in the understory to promote sugar maple recruitment. 75 From a slightly longer-term view, stand-level median sapling ages are currently 66 years old; these stems would be expected to reach economic maturity in 40 years. Although possible, this scenario would require full release of stems, and their quality may be compromised due to extended suppression. Furthermore, saplings are currently understocked (Figure 3.1), so their release would likely result in perpetual understocking of larger size classes. Stands with most saplings dominated by old, low-vigor saplings may benefit from eliminating or reducing current sapling stock and attempting to promote high-light growth environments to release younger and presumably more vigorous advanced regeneration. From a regional (Michigan) perspective, alternative silvicultural regimes may need to be implemented imminently across broad areas to diversify sugar maple age structure (Appendix D, Figure D.3); this may be necessary to minimize disruption to economic supply of sugar maple timber products, as many of these stands will be approaching economic maturity simultaneously. 3.5.4 Caveats and future directions An important caveat for our method of identifying age cohorts was that it was biased to minimize age cohort width and maximize number of age cohorts to rigorously test our hypothesis that most stands would have few age cohorts. Some identified cohorts close in age may truly represent the same recruitment event and therefore age cohort, which would further support our notion that recruitment or establishment events are likely to be diffuse in age range. Furthermore, some potential low-density, old-age (e.g., > 200 years) age cohorts are missed by our methodology by employing a single threshold across all diameter classes, although stands with older trees overall tended to have more empty diameter bins (0 stems ha-1) and therefore a lower threshold relative to total sugar maple stem density. 76 Future analyses, both with our samples and in other stands, can contribute significantly to our understanding of sugar maple age structure in STS managed NHF. Additional analyses with more refined measures (e.g., soil measurements for site quality) may shed additional insights on the influence of deer browsing and nutrient and water availability in determining sugar maple age structure. Further analysis with a subset of this data to characterize ring width over time would enable a look at past suppression and recruitment dynamics. This could yield insight into how recruitment patterns have changed over time, such as average age and diameter of stems at canopy recruitment. Additionally, more intensive, directed sampling of stems 5 – 10 cm as well as 0 – 5 cm could yield more refined insight into the age – diameter relationships in smaller diameter classes. We assumed sampled saplings represented all saplings within the 5 – 10 cm class, but additional sampling intensity may reveal greater age diversity; given the highly suppressed nature of saplings, 5 cm DBH vs. 10 cm DBH stems may vary widely. Additional visual characterizations of these small stems may also reveal traits that may serve as a better proxy for age than diameter, height, or live crown ratio, such as crown shape (Bartholomew, 2018). 77 CHAPTER 4 PATTERNS AMONG SPECIES FROM SEEDLING TO CANOPY PORTEND FUTURE COMPOSITIONAL SHIFTS AND LOW RESILIENCE IN MANAGED NORTHERN HARDWOOD FORESTS 4.1 Abstract Declining biodiversity is pervasive globally, challenging sustainability and resilience of ecosystems facing multiple environmental stresses. Ecologically and economically valuable northern hardwood forests of the Great Lakes region have experienced declines in canopy tree diversity since European colonization. Current tree regeneration trends (i.e., restricted to shade tolerant, deer browsing resilient/resistant species) in managed stands suggest low tree-diversity patterns among mature/canopy trees will continue and may intensify in the future; however, studies characterizing current diversity and projecting further trends at regional (106 ha) scales are lacking. For 141 managed northern hardwood forests across northern Michigan, USA, we analyzed patterns and possible drivers of species diversity and relative abundance by species from seedlings to canopy trees. Overall, the estimated asymptotic stand-level richness for seedlings (< 137 cm tall) and saplings (> 137 cm tall and <10 cm diameter) averages 6 and 7 species, respectively, and with relatively low evenness (effective number of species, Hill q1, < 3 species). Stands with lower species diversity of canopy/mature stems (> 10 cm DBH) tend to have lower seedling and sapling diversity. Models of individual species abundance indicate that high total canopy/mature basal area (i.e., shade) and low conspecific mature stem density (lack of seed sources) consistently contributed to low seedling/sapling density. Among the 17 most common species, seven have lower relative importance (RI) as seedlings than mature trees (e.g., Tsuga canadensis, Betula alleghaniensis, Populus grandidentata, Tilia americana) suggesting 78 seedling establishment (germination substrate and/or early survival) limitations; two (Acer saccharum, Quercus rubra) had lowest RI as saplings suggesting seedling to sapling recruitment bottlenecks (e.g., deer browsing); six (e.g., Fagus grandifolia, Ostrya virginiana, Pinus strobus) had highest RI as saplings indicating low seedling to sapling recruitment bottlenecks; two (Prunus spp., Acer rubrum) had highest RI as seedlings with moderate declines in sapling and canopy/mature classes indicating modest sapling and canopy/mature recruitment bottlenecks; and Fraxinus americana had highest RI as seedlings and lowest for canopy/mature trees, which is consistent with its recent widespread canopy mortality from emerald ash borer. Observed regeneration diversity patterns suggest future shifts in mature canopy stratum composition and highlight a need for changes in management aimed at increasing tree species diversity and desirable species mixes. 4.2 Introduction Biodiversity contributes to ecosystem resilience, or the capacity of a system to recover following disturbance. Resilience via diversity will be necessary for ecosystem functioning in the face of increasingly frequent and diverse human-mediated disturbances including forest fire activity (Flannigan et al., 2000), invasive pests and pathogens (Dukes et al., 2009; Ramsfield et al., 2016), and frequency and severity of drought (Allen et al., 2010), plus interactions among factors (Dale et al., 2000) and with climate change (Dale et al., 2001). Unfortunately, for North American temperate forest ecosystems, accelerating multiple threats to ecosystem functioning comes at a time of declining tree species regeneration diversity (Ramirez et al., 2018; Miller & McGill, 2019; Vickers et al., 2019), in forests already documented as having declining canopy diversity (Schulte et al., 2007); these patterns that may portend future stasis or further declines in canopy biodiversity and resilience. Characterizing tree regeneration species diversity broadly, 79 geographically, will be important for forecasting future forest resilience globally under increasingly volatile environments. Covering over 50 million hectares of northeastern United States (Oswalt et al., 2014), northern hardwood forests (NHF) are an important forest ecosystem currently dominated by ecologically and economically valuable sugar maple (Acer saccharum Marsh) (Linehan & Jacobson, 2005; Duval et al., 2014). Compared to pre-European colonization, current canopy tree species richness and structural complexity are lower overall (Schulte et al., 2007) and may be lower in managed than unmanaged stands (Crow et al., 2002; Neuendorff et al., 2007; Powers & Nagel, 2009). Low diversity and density in the regeneration layer of key species, such as sugar maple, are particularly concerning (Leak, 2006; Neuendorff et al., 2007; Powers and Nagel, 2009; Matonis et al., 2011). However, few studies have analyzed stand-level tree regeneration diversity at regional (106 ha) scales, which may capture driving factors that operate at that scale. Addressing challenges to securing diverse tree regeneration is becoming increasingly important, given predicted declines in ecosystem functioning in the Great Lakes region from climate change (Rogers et al., 2017) and continued invasive pest and pathogen establishment and spread (Liebhold et al., 2013; D. Chapman et al., 2017; Panzavolta et al., 2021). Although diversity alone can contribute to resilience in the face of multiple potential disturbances, species composition also matters, and forests may be undergoing compositional shifts that change their capability to provide multiple ecosystem goods and services that are resilient to diverse disturbances (Walters et al., 2022). Many studies indicating declining tree regeneration density for key species also show, or suggest, different composition for regeneration than canopy trees (e.g., Neuendorff et al., 2007; Matonis et al., 2011; Elenitsky et al., 2020; Hupperts et al., 2020), indicating possible compositional shifts are underway. However, there are 80 few explicit examinations of composition in seedling, sapling and canopy classes and factors potentially underlying differences across size class. Quantifying current northern hardwood tree regeneration composition and diversity and associated driving factors is key for forecasting potential future canopy composition and diversity and for development and strategic implementation of adaptive forest management. Several potential driving factors likely influence tree regeneration composition and diversity in managed northern hardwood forests. Abiotic drivers of nutrient and water availability, which vary spatially, influence tree species distribution, assemblages, and therefore diversity. For the Great Lakes region, post-glacial landforms, which vary in nutrient availability and soil water holding capacity (site quality), help drive forest type distribution, from xeric jack pine-black oak dominated coarse textured outwash soils, to mesic northern hardwood dominated finer textured moraine soils (Zak et al., 1986, 1989; Baribault et al., 2010). Interacting with site quality is a Great Lakes-climate driven (i.e., lake-effect) three-fold gradient of annual snowfall, further influencing forest distribution and species distributions, particularly for sugar maple (Henne et al., 2007). Over the narrower range of higher site qualities supporting northern hardwood forests, tree regeneration composition varies (Elenitsky et al., 2020), but the potential for site to influence future composition of all strata could be overridden by other factors, such as deer browsing or management (Matonis et al., 2011; Bannon et al., 2015). Compositional data supporting habitat classification development for the region (i.e., frequency of occurrence of saplings, Burger and Kotar 2003), suggest that lower site quality sites supporting NH tend to have more species; however, explicit examinations of tree diversity in relation to site quality are lacking. 81 Interacting with the abiotic landscape, biotic drivers contribute to patterns of tree diversity in NHF of the Great Lakes region. Invasive pests and pathogen epidemics have decreased the likelihood of canopy ascendance/maturity for American beech (Fagus grandifolia), white ash (Fraxinus americana) and American elm (Ulmus americana) (Parker & Leopold, 1983; Forrester et al., 2003; Nuckolls et al., 2009; Klooster et al., 2014). Impacting a wider range of tree species, invasive non-native earthworms alter nutrient cycling, retard decomposition, and reduce understory plant abundance (Bohlen et al., 2004; Holdsworth et al., 2007, 2008; Resner et al., 2015), potentially influencing tree regeneration composition (Hale et al., 2006). Abundant white-tailed deer (Odocoileus virginianus) populations, which have increased drastically since the early 1900s (MDNR, 2016), can shift species composition through selective browsing and depress overall regeneration density (Côté et al., 2004; Bradshaw & Waller, 2016), including canopy dominant sugar maple (Curtis and Rushmore, 1958; Horsley et al., 2003; Kain et al., 2011; Matonis et al. 2011; White, 2012, Henry et al. 2021). Deer browsing interacts with stand level variables to influence tree regeneration density and composition, including site quality (Randall & Walters, 2011) and light availability/gap size (Sage et al., 2003; Walters et al., 2016). Since deer browsing is vertically limited to a maximum of ~2 m (Walters et al., 2020a), regeneration which has vertically escaped the deer browsing zone may be particularly indicative of future canopy structure and composition. Capturing variability in biotic drivers, such as invasive pests and pathogens or local deer populations, is key to fully characterizing trends in tree regeneration diversity and distribution, which is necessary for identifying particularly low-diversity regions in need of alternative management intervention. In addition to biotic and abiotic drivers, past forest management may have played both direct and indirect roles in patterns of tree regeneration diversity and density. The single-tree 82 selection (STS) silvicultural system was introduced in the late 1950’s (Arbogast 1957) and has since dominated northern hardwood forest management (Kern et al., 2014). Under this system, periodic partial harvests every 10-20 years remove select trees, resulting in small, dispersed harvest gaps. These gaps theoretically promote the regeneration of economically desirable shade- tolerant species, such as sugar maple. Under this system, tree diversity may decline via both intentional removal of competing stems of less economically valuable species and decreased regeneration and canopy recruitment of less shade tolerant species in small STS gaps (Niese & Strong, 1992; Crow et al., 2002). Local (stand) scale, species-specific declines or elimination of mature/canopy tree seed sources can negatively impact regeneration rates for many NHF species (Willis et al, 2016), particularly those with shorter distance dispersal potential (McEuen & Curran, 2004). Small seeded species, such as yellow birch and eastern hemlock, may face additional seedling substrate establishment constraints via low coarse wood coverage resulting from periodic harvests capturing potential large tree mortality and from winter logging and use of lower impact harvesting equipment resulting in minimal exposure of mineral soil substrate via disruption of the forest floor (Marx & Walters, 2008; Bolton & D’Amato, 2011; Willis et al., 2016). Given widespread implementation of STS management in NHF, characterizing tree diversity regeneration outcomes for these forests is key for future management success as well as landscape-level forest health and functioning. Despite potential failure of STS to promote diverse forests, management can play a role in improving future diversity. Climate change forest management will involve efforts to adapt, mitigate, or migrate forests (Aitken et al., 2008), and managers will need to be increasingly flexible, employing a variety of strategies to adapt to future conditions (Millar et al., 2007). For northern hardwood forests, alternative, more intense management harvesting regimes are being 83 explored (Hupperts et al., 2020; Walters et al., 2020; Rogers et al., 2021). These include testing the notion that larger canopy harvest gaps promote more species diverse tree regeneration by promoting admixtures of shade-intolerant tree species (e.g., Niese and Strong, 1992; Walters et al., 2016). However, some studies suggest more intense harvests may not increase regeneration diversity in the long-term (Kern et al., 2013; Knapp et al., 2021; but see Niese and Strong, 1992). Accurately characterizing current patterns and associated driving factors is necessary to guide regional implementation of future management of NHF, whether by continued use of STS or by alternative silvicultural strategies. In this study we address knowledge gaps in understanding the patterns and potential causes of current and potential future diversity and composition of NHF. To do this we quantified patterns of individual species regeneration abundance and diversity, plus associations with likely stand and landscape drivers. We also characterized variation in species diversity and abundance as a function of size class and interpret these patterns in terms of potential drivers of pattern and implications for compositional shifts in the future. We address these goals by analyzing a dataset of 141 selection-managed northern hardwood forests, distributed throughout northern Michigan, which capture a variety of potential landscape and stand level drivers. Given current management regimes and other abiotic and biotic factors, we predict: 1) Tree diversity will be low particularly on high quality sites and in all size strata, but especially low in seedling and sapling layers due to filtering effects of substrate, resource limitation, and selective deer browsing pressure. 2) As driven by limiting factors listed in (1) and others, we predict compositions in seedling, sapling, and canopy strata will differ, portending possible future shifts in canopy composition. 84 2a) Canopy composition (conspecific basal area of stems > 10 cm DBH) will associate with regeneration density and diversity because of seed source limitation, particularly for large-seeded species 2b) Small seeded species will face seedling establishment limitations, negatively associating with hardwood litter coverage, and having greater relative representation in the overstory than in the understory 2c) Shade intolerant and mid-tolerant species will associate negatively with increasing basal area and will have greater representation in the canopy/mature strata versus the understory 2d) Highly palatable tree species, like sugar maple, will negatively associate with deer use for the sapling class and will have greater representation in the canopy or seedling layers versus the sapling layer We interpret these results in the context of implications for future tree diversity patterns and with explicit consideration to management implications. 4.3 Methods 4.3.1 Study area and data collection Our data are from 141 managed northern hardwood forest stands distributed throughout northern Michigan which are part of a larger experiment analyzing harvest methods and regeneration outcomes (Walters et al. 2020b). Here, we focus on pre-harvest regeneration. These stands were considered ready or near ready for partial harvest by single-tree selection criteria (i.e., > 23 m2 ha-1 BA and well stocked in sawtimber classes, Arbogast, 1957) at the time of sampling and were likely subject to 3-4 partial/selections harvests prior to our sampling them. 85 Additional vegetation sampling details can be found in Henry et al. 2021); we recount pertinent details here. Most stands were state owned (n=119) and managed by the Michigan Department of Natural Resources (MDNR), with the remaining 22 stands owned and managed by private forest industry companies (Hancock Timber Resource Group, The Rohatyn Group). Stands were generally on upland sites, mesic to wet mesic, with fertile medium-textured upland soils (Dickmann & Leefers, 2016) and dominated by sugar maple (vast majority), red maple or basswood (see Table 4.1 for scientific names and authority). We established a systematic grid of 25 survey points within a 12-ha square within each stand and collected species-level data in the following plots at each survey point. We counted all tree seedlings shorter than 50 cm and ocularly estimated percent cover of shrubs and groundcover composition (e.g., hardwood litter, bare mineral soil, downed woody material) within 2 quadrats, each 1 m2 (50 m2 total per stand). Within a 2 m radius circular plot (314 m2 total per stand) we tallied tree regeneration for stems 50 to 137 cm tall and tallied and measured DBH for stems > 137 cm tall and < 5 cm DBH. Within a 6 m radius circular plot (2,827 m2 total per stand), we tallied and measured DBH for all stems > 5 cm DBH. Diameters were converted to basal areas (stem cross sectional areas at breast height), summed, and expressed on a m2 per ha basis, with values (BA) representing the density of trees > 10 cm DBH (i.e., superior to our largest sapling class). Using herbaceous plant assemblages and key indicator species, we assigned each stand a dominant habitat type, which is an index of site quality (Burger & Kotar, 2003). We assessed winter deer use with fecal pellet transect surveys in spring of 2017 at a subset of 50 sites, surveying 628 m of transect (6 m wide) subdivided into 103 segments at each site. We modeled estimated winter deer use (percent transect segments occupied by deer pellets) 86 for all stands in 2017 with a spatial model in a Bayesian framework; the model also incorporated well-established climate and landscape factors known to influence local deer populations, plus deer pellet transect surveys conducted in 2019 at 139 sites following experimental harvest (Appendix A). 87 Common Mean stems ha-1 (Max stems ha-1) | n Scientific name CODE name Seedlings Saplings Canopy n Abies balsamea [L] ** BF Balsam fir 817 (11,295) | 47 154 (936) | 64 41 (279) | 62 77 Mill. Acer rubrum L. ** RM Red maple 17,823 (139,733) | 118 641 (5,673) | 71 86 (460) | 84 119 Acer saccharum Marsh. ** SM Sugar maple 71,797 (448,752) | 140 1,217 (6,943) | 140 282 (531) | 141 139 Betula alleghaniensis ** YB Yellow birch 609 (3,664) | 47 81 (725) | 39 15 (113) | 66 85 Britton Betula papyrifera * PB Paper birch 350 (1,055) | 9 10 (32) | 7 10 (28) | 20 27 Marsh. Fagus grandifolia ** AB American 4,090 (49,171) | 89 1,671 (5,903) | 96 41 (195) | 87 98 Ehrh. beech Fraxinus americana L. ** WA White ash 7,885 (55,453) | 71 451 (5,188) | 43 12 (67) | 31 72 Fraxinus nigra Marsh. BA Black ash 64 (64) | 1 39 (39) | 1 4 (4) | 5 5 Fraxinus pennsylvanica GA Green ash 1,483 (5,600) | 5 35 (64) | 3 9 (18) | 5 8 Marsh. Ostrya virginiana Mill. ** IW Ironwood 2,744 (14,566) | 108 955 (4,545) | 111 17 (74) | 77 116 Picea spp. * SP Spruce species 244 (2,159) | 22 25 (223) | 35 9 (39) | 32 3 Table 4.1. Summary statistics for tree species regeneration at the stand level when present in a given size class, including mean, maximum (stems ha-1), and number of stand occurrences (n, out of 141) by size class, plus overall number of site occurrences n (across all size classes). Summary statistics were calculated for the subset of stands which that species was present on, in a given size class. Seedlings are 0 – 137 cm tall, small saplings are > 137 cm tall and < 10 cm DBH, and canopy trees are > 10 cm DBH. Species labeled with ** were modeled for individual abundance and relative abundance by size class; species labeled with * were modeled only for relative abundance by size class. Species without any * were included only for general species diversity calculations and stand summaries. 87 Figure 4.1[cont’d] Common Mean stems ha-1 (Max stems ha-1) | n Scientific name CODE name Seedlings Saplings Canopy n Pinus resinosa Ait. RP Red pine NA (0) | 0 4 (4) | 1 51 (85) | 2 53 Pinus strobus L. * WP White pine 279 (2,059) | 22 77 (700) | 17 9 (21) | 14 30 Populus balsamifera L. BP Balsam poplar NA (0) | 0 187 (187) | 1 NA (0) | 0 1 Populus grandidentata * BTA Big tooth aspen 459 (1,495) | 5 73 (127) | 7 12 (67) | 13 20 Michx. Populus tremuloides * QA Quaking aspen 342 (833) | 9 98 (421) | 11 11 (57) | 23 31 Michx. Prunus pennsylvanica * PC Pin cherry 497 (1,632) | 6 70 (396) | 10 6 (11) | 3 16 L.f. Prunus spp. ** CH Cherry species 3,283 (31,632) | 112 126 (863) | 73 32 (191) | 59 11 8 Quercus rubra L. ** RO Red oak 1,166 (3,382) | 35 138 (605) | 13 16 (64) | 26 37 Sorbus americana MA Mountain ash 505 (1,600) | 6 4 (4) | 2 11 (11) | 1 7 Marsh. Thuja occidentalis L. EWC Eastern white cedar 351 (805) | 4 115 (226) | 2 18 (53) | 8 10 Tilia americana L. ** BW Basswood 884 (3,832) | 65 89 (732) | 31 49 (279) | 78 87 Tsuga canadensis [L.] * EH Eastern hemlock 253 (864) | 13 44 (138) | 13 16 (88) | 29 35 Carr. Ulmus americana L. ** AE American elm 441 (2,400) | 31 47 (389) | 20 7 (21) | 20 47 Ulmus rubra Muhl. SE Slippery elm 216 (432) | 4 33 (35) | 3 4 (4) | 3 4 88 4.3.2 Species diversity analyses All analyses and modeling were conducted in R (R Core Team, 2018). To characterize species diversity, we chose to calculate Hill numbers, a unified set of diversity measures parameterized by a diversity order q, which correspond to common diversity metrics (species richness, q=0; exponential of Shannon value, q=1; the inverse of Simpson diversity, q=3). All Hill numbers are interpretable as number of species; for example, q = 0 is the number of species present, and q = 1 is the effective number of common species, or the number of species, equally abundant, which would result in the same Shannon value (Chao et al., 2014). Given low species richness of our stands, we chose to focus on Hill numbers q = 0 (species richness, SR) and the q = 1 (effective number of common species, CS) (Hsieh et al., 2016). We were limited in having sampled a fixed area for all size classes, although we did sample larger size classes in larger sample plots to attempt to achieve similar sampling intensity by size class. Despite this, there could still be imbalances in sampling intensity (by number of individuals sampled) among size classes and sites due to variable stem density. Because of these considerations, we used the iNEXT package (Hsieh et al., 2016), which conducts rarefaction and extrapolation as outlined in Chao et al. (2014) to estimate asymptotic Hill numbers. Rarefaction accounts for sampling completeness and involves generating a rarefaction curve, which shows number of species as a function of sampling intensity (number of individuals or number of sampling units). iNEXT conducts rarefaction and then extrapolates the curve to estimate the value at which the curve asymptotes, representing the theoretical species richness as sampling effort achieves completeness. We analyzed data as abundance values (total number of individuals sampled, per species) rather than incidence values (total number of plot occurrences 89 per species), given the relatively small number of plots (25). We used the iNEXT function to calculate SR and CS for all three size classes. We conducted generalized linear models for both diversity metrics, by size class, as a function of landscape and site-level factors, which we standardized for modeling: winter deer use (deer), average stand-level basal area (BA, m2 ha-1), standard deviation of stand-level basal area (BASD), canopy/mature (> 10 cm DBH) asymptotic species richness (SR10) and effective common species (CS10), plus groundcover percentages of hardwood litter (HWL), bare mineral soil (BMS), and coarse woody debris (CWD) and shrub coverage (shrub) for the seedling class models (Table 2.2). While we focus on the two most important groundcover classes (HWL and BMS), there were a variety of other groundcover classes recorded, such that BMS and HWL do not sum to 1. We calculated variance inflation factors (VIF), a measure of multicollinearity among predictors in linear regression models (Marquaridt, 1970), using check_collinearity in the performance package (Lüdecke; Makowski; Waggoner; Patil, 2020); at a threshold of 5, below which predictors are generally considered to have low collinearity (James et al., 2013), none warranted removal. In a Bayesian framework, we used a gamma distribution to model estimated asymptotic SR and CS (both of which are continuous positive). We ran models using the rjags package (Plummer, 2018), which interfaces JAGS (Just Another Gibbs Sampler) software (Plummer, 2003) with R. We ran models for 20,000 iterations after a 5,000-step adaptation period with three chains, which, based on our visualizing trace plots and calculating Gelman’s Diagonal, was sufficient for model convergence. We used the coda package for model convergence diagnostics (Plummer et al., 2006). 90 Predictor Acronym Mean Min Max Deer use (% segments occupied) DEER 19 2 60 Basal area (m2 ha-1), stems > 10 BA 28 20 46 cm DBH Basal area standard deviation BASD 14 7 23 -1 Conspecific density (stems ha ) SeedS See Table 4.1 Canopy of stems > 10 cm DBH † Asymptotic species richness, SR10 7 1 20 stems > 10 cm DBH Hill q=1 (effective common CS10 3 1 7 species) for stems > 10 cm DBH Shrubs (% cover) * SHRUB 2 0 26 Bare mineral soil (% cover) * BMS 1 0 56 Hardwood litter (% cover) * HWL 91 27 97 Coarse woody debris (% cover) * CWD 3 1 7 Table 4.2. Summary statistics for predictor variables included seedling and sapling models for diversity and individual species distribution. Variables include estimated deer use (DEER), average basal area (m2ha-1, BA), standard deviation of basal area (BASD), seed source as estimated by the density of conspecific stems > 10 cm DBH (SeedS), estimated asymptotic species richness of stems > 10 cm DBH, estimated effective number of common species (Hill q=1) of stems > 10 cm DBH, % coverage of shrubs (shrub), % coverage of bare mineral soil (BMS), hardwood litter (HWL), and coarse woody debris (CWD). † SeedS was only included in the individual species models. * SHRUB, BMS, HWL, and CWD were only used as predictors for the seedling models. We evaluated associations with site quality through a separate analysis to preserve degrees of freedom. To determine whether ordinal site quality diversity measures were comparable between the three broad geographic regions of our study, we compared generalized linear models of diversity predicted by site quality versus a site quality – region interaction, using gamma distributions. We compared the simple model versus the interaction model using the Akaike information criterion (AIC), which is an estimator of prediction error (Akaike, 1974). The interaction model had a lower AIC for all three size classes and both diversity metrics. Since we were primarily interested in the effect of site quality, we used the emmeans function to 91 estimate marginal means and the pairs function to conduct pairwise comparisons across site quality for each region separately (Lenth, 2020). 4.3.3 Generalized linear modeling for species abundance We included 11 species present in > 10% of stands (> 14 stands) for each of the three size classes. We included the same predictors as in the species diversity models (Table 4.2), plus an additional predictor to measure density of conspecific stems > 10 cm DBH (SeedS), a proxy for seed source (McEuen & Curran, 2004). We did not include site quality as a predictor because our chosen index of site quality and tree composition are confounded and because tree composition is included in the model in the SeedS predictor (Burger & Kotar, 2003). For each species, we subset the data to stands which that species was present on, in any size class. We compared two candidate models: a negative binomial generalized linear model and a zero inflated negative binomial generalized linear model. We checked for multicollinearity, simulated models in a Bayesian framework, and checked for model convergence using the same procedures as described for species diversity metric models. To select a model, we compared Deviance Information Criteria (DIC) values, which incorporates model fit and complexity (Spiegelhalter et al., 2002). DIC comparison led to model selection of a negative binomial model for seedlings and saplings of AB, CH, IW, RM, and SM, plus BW seedlings. BW saplings, plus AE, BF, RO, WA, and YB were modeled with a zero-inflated negative binomial model following DIC comparison (Appendix E, Table E.1). We reported selected model results at p > 0.10 given relatively low number of stands, coarse aggregation of data, and our intention to highlight trends. 4.3.4 Relative abundance as a function of size class For each stand, we calculated relative abundance for each present species by size class (stems of a given species divided by total stems). We limited our relative abundance by size class 92 analyses to the 18 tree species present, in any size class, on at least 10% of stands (> 14 stands). In this case, relative abundance is useful because all three size classes are portrayed on the same scale, independent of stem density, which varies widely between seedlings and canopy/mature trees; it also has the benefit of reflecting potential shifts in composition among classes as it places species-specific density in the context of competitors. For each stand, we extracted the rank order of species abundance, by size class, for each species present. For example, in a hypothetical stand, sugar maple has the highest relative abundance in canopy/mature (referred to as canopy for brevity) stems (0.7), followed by seedlings (0.6) and saplings (0.4), making the ranking “canopy – seedling – sapling”; in the same stand, hemlock had a relative abundance of 0.01 in the canopy, and 0 for seedlings and saplings, making the ranking just “canopy”. We then tallied the frequency of each relative abundance ranking among stands, for each species. For example, sugar maple had the ranking “canopy – seedling – sapling” in X stands, “seedling – canopy – sapling” in Y stands, and so on. We only included sites in which a species was present in at least one size class for this analysis. This resulted in n < 141 for all species except sugar maple, which was present on all sites. We used a Kruskal-Wallis rank sum test to determine if relative abundance varied by size class, and, if significant, pairwise Wilcoxon tests for paired (within stand) values with a Bonferroni correction (n=48) to determine individual differences. To summarize differences in relative abundance among size class by species, we labeled seedling- sapling-canopy distributions with four shape categories: left skewed, right skewed, positive parabola, and negative parabola. 93 4.4 Results 4.4.1 Overview of stand characteristics BA of trees > 10 cm DBH averaged 28 m2 ha-1, and BASD averaged 14 m2 ha-1 (Table 4.2). SHRUB, BMS, and CWD were all low percentages (< 3 %), compared to HWL (91%). DEER was variable, ranging from 2% to 60% pellet occupancy of transect segments. We documented 25 tree species in our survey. Total stem density (stems ha-1) in stands ranged from 9,605 – 460,112 (mean 99,033) for seedlings, 357 – 12,771 (mean 3,797) for saplings, and 219 – 778 (mean 449) for canopy/mature stems (>10 cm DBH). Nine species were sampled in any size class on at least half of the sites, including sugar maple which was sampled on all stands (Table 4.1). Summing across all 141 stands, nine species had an average relative abundance of at least 1% in at least one size class (Figure 4.1). Although sugar maple had the highest average relative abundance for seedlings, saplings, and canopy/mature trees, relative abundance of ironwood and American beech were nearly equivalent to sugar maple for the sapling class (Figure 4.1). 94 Figure 4.1. Landscape-level (i.e., across 141 stands in Michigan) average relative abundance, by size class, for species which achieved at least 1% relative abundance in any size class; species include sugar maple (SM), American beech (AB), ironwood (IW), red maple (RM), balsam fir (BF), white ash (WA), cherry species (CH), basswood (BW), and yellow birch (YB). Scientific nomenclature and all species included in “Other” can be found in Table 4.1. Seedlings are 0 – 137 cm tall, saplings are > 137 cm tall and < 10 cm DBH, and canopy/mature are > 10 cm DBH. Among stands and size classes the highest estimated asymptotic species richness was 20 and the lowest was 1 (Figure 4.2). Among size classes, average species richness was highest for seedlings (7.6), followed by canopy/mature stems (7.4) and saplings (6.4), though only seedlings and saplings differed from one another (p < 0.05) (Figure 4.2). Compared to species richness, effective common species was generally lower and averaged 3.2 for saplings, 2.8 for canopy/mature trees, and 2.4 for seedlings, with all size categories different from one another (p < 0.05) (Figure 4.2). 95 Figure 4.2. Stand-level estimated asymptotic species richness (SR) and effective number of common species (CS) for seedlings (0 –137 cm tall), saplings (> 137 cm tall and < 10 cm DBH), and canopy/mature trees (abbreviated canopy for brevity, > 10 cm DBH). Asymptotic estimates account for sampling effort and are the estimated true number of species present in the stand. Within each panel, box plots labeled with the same letter are not significantly different. 4.4.2 Species diversity measures predicted by stand and landscape factors There were few significant associations of species richness and diversity with factors we assessed (Table 4.3). Seedling diversity (SR, CS) was generally positively associated with canopy/mature diversity (CS10), and saplings SR was negatively associated with BA. For substrate, seedlings SR declined with BMS, while CS decreased with HWL. Diversity measures were generally flat or declining from low to high site quality for all size classes, with declines only significant for CS for seedling and canopy/mature strata in the Eastern Upper Peninsula (Figure 4.3). 96 Size SR CS Variable class 50% 95% CI 50% 95% CI Intercept 2.3106 (1.3203 - 3.348) 1.856 (0.7659 - 2.986) Shrub 0.0069 (-0.0086 - 0.0234) -0.0102 (-0.0274 - 0.0083) Deer 0.0000 (-0.0054 - 0.0055) -0.0004 (-0.0067 - 0.0062) BA 0.0027 (-0.0121 - 0.0177) 0.0004 (-0.0158 - 0.0159) BASD 0.0089 (-0.0124 - 0.0297) 0.0095 (-0.0145 - 0.0332) Seedlings BMS -0.0204 (-0.0352 - -0.0047) -0.0079 (-0.0243 - 0.0099) HWL -0.0079 (-0.0184 - 0.002) -0.0121 (-0.0225 - -0.0019) CWD -0.0158 (-0.0672 - 0.0357) -0.0355 (-0.0971 - 0.0261) SR10 0.0251 (0.0062 - 0.0446) -0.0077 (-0.0304 - 0.0158) CS10 0.0498 (-0.0077 - 0.1075) 0.0986 (0.0287– 0.1693) Beta 9.5024 (7.4123 - 11.9492) - Shape - 7.0027 (5.4912 - 8.7879) Intercept 1.8002 (1.3858 - 2.2575) 0.9976 (0.6065 - 1.4053) Deer -0.0035 (-0.0089 - 0.0019) -0.0007 (-0.0057 - 0.0044) BA -0.0156 (-0.03 - -0.0016) -0.0093 (-0.0224 - 0.0038) Saplings BASD 0.0052 (-0.0171 - 0.0271) 0.0034 (-0.016 - 0.0235) SR10 0.0125 (-0.0076 - 0.033) 0.0082 (-0.0096 - 0.0268) CS10 0.1308 (0.07 - 0.1921) 0.1111 (0.0558 - 0.1669) Beta 8.4001 (6.6 - 10.5038) - Shape - 9.9588 (7.7741 - 12.5021) Table 4.3. Parameter estimations for models of species richness (SR) and diversity (CS) for seedlings (0 – 137 cm tall), saplings (> 137 cm tall and < 10 cm DBH), and canopy/mature trees (> 10 cm DBH), as predicted by total canopy/mature basal area, standard deviation of canopy/mature basal area, and deer use, plus shrub coverage and percent hardwood litter coverage (HWL), coarse woody debris (CWD), and bare mineral soil (BMS) for the seedling models. Model fit parameters (Intercept, Beta and Shape) are also provided: Beta characterizes dispersion in a negative binomial distribution (SR model only) and shape influences the shape of the gamma distribution (CS model only). Significant associations (95% CI) are bolded. 97 Figure 4.3. Histograms for asymptotic species richness (SR) and number of common species, CS, by site quality and region. Site quality 1 = lowest quality and 4 = highest (Burger and Kotar 2003). Letters denote significant differences (pairwise comparison of estimated marginal means) among site qualities within individual panels; bars with the same letters are not significantly different. 98 4.4.3 Individual tree species density predicted by stand and landscape factors In contrast to community-level species richness/diversity metrics, seedling and sapling densities of individual species had numerous significant stand-level associations with key drivers, including Deer, BA, and SeedS (Figure 4.4). Balsam fir seedling densities were positively associated with BA, while saplings of red maple, sugar maple, and cherry species negatively associated with BA. Maple species had higher densities in stands with greater BA SD, while balsam fir had negative associations. SeedS had, by far, the most consistent associations among factors potentially driving species specific densities, with positive associations for all species in both seedling and sapling classes. Biotic factors were significant in a few cases, including sugar maple sapling density negatively associating with DEER and positively with SHRUB. Lastly, substrate had few significant associations with seedlings (the only size class tested), but they lacked generalization, again perhaps, due to limited variation and low coverage of CWD and BMS substrates. 99 Figure 4.4. Parameter estimates (median with 0.05 – 0.95 confidence interval) for seedling (0 – 137 cm tall) and sapling (> 137 cm tall and < 10 cm DBH) densities with predictor variables by tree species: average stand canopy/mature basal area (BA), standard deviation of canopy/mature basal area (BASD), deer use (Deer), seed source as estimated by conspecific canopy/mature tree density (SeedS), shrub coverage (Shrub), percent bare mineral soil (BMS), percent hardwood litter (HWL), and percent coarse woody debris (CWD). CWD, HWL, BMS, and Shrub were only included in the seedling model. Median values of significant (p > 0.10) parameters are diamonds, outlined in black. Seedlings and saplings of AB, CH, IW, RM, and SM, plus BW seedlings were modeled with a negative binomial model; all others were modeled with a zero-inflated model. 100 4.4.4 Species relative abundance across size classes Of the 18 species analyzed, relative abundance varied significantly by size class for 16 species (Figure 4.5). Size class distributions were left-skewed (relative abundance increasing as a function of size) for seven species (e.g., eastern hemlock), a negative parabola shape (saplings having highest relative abundance) for six species (e.g., American beech), positive parabola (saplings least abundant class) for two species (e.g., sugar maple), and right-skewed (relative abundance declining with size class) for three species (e.g., white ash). 101 Figure 4.5. Relative abundance by species relative to stems of all species of the same size class for seedlings (0 –137 cm tall), saplings (> 137 tall and < 10 cm DBH), and canopy/mature (abbreviated canopy in figure, > 10 cm DBH) strata. Within a species plot, different letters indicate different significances in relative abundance between size classes (p-value < 0.05). Size class was not a significant predictor for American elm and white pine, so we did not calculate pairwise comparisons. Histogram color corresponds to their identified shape across size classes (legend on right). 102 Summaries of relative abundance rankings among size classes for each stand (Table 4.4). are generally consistent with patterns suggested by the average, landscape-level patterns of relative abundance (Figure 4.5); however, stand by stand variation reveals additional details important in terms of species ecologies and management (see Discussion). For example, left skewed species are defined by canopy/mature trees as the dominant relative abundance category (Figure 4.5), and this pattern is generally supported by the frequency of relative abundance patterns at the level of individual stands (for 1st and 2nd most common patterns, canopy trees ranked first, Table 4.4). However, for individual stands, it is striking that for four of the seven species (yellow birch, paper birch, big tooth aspen and eastern hemlock) the most common stand/site level relative abundance pattern was to have only canopy/mature trees in the sample (i.e., no saplings or seedlings). Similar but oppositely for right-skewed species for two of the three species (white ash, cherry), the most common relative abundance pattern at the stand level was for samples to contain only seedlings. There were also important differences within shape groups. For example, over the landscape, red oak and sugar maple have lower relative abundance of saplings than seedlings and canopy/mature trees (positive parabola), but saplings are still present (relative abundance = 0.32 for SM and 0.01 for RO). Among stands, the most common relative abundance pattern reflects this overall landscape pattern (relative abundance of seedling > canopy > saplings). However, for red oak, the most common pattern was to have only seedlings. Although stand-level patterns (particularly size-class absences for a given species) are influenced by stem density and sampling completeness, they have clear broad implications for species ecologies and management (see Discussion). 103 Seedling Sapling Canopy 1st ranking 2nd ranking Shape Species n % % % Order n Order n BF 77 0.09 0.32 0.58 Can-Sap-Seed 21 Can-Sap 15 YB 85 0.16 0.14 0.69 Can 23 Can-Seed 14 PB 27 0.19 0.07 0.74 Can 13 Seed 5 Left SP 53 0.11 0.4 0.49 Sap 14 Can 9 skewed BTA 20 0.15 0.3 0.55 Can 8 Sap 4 BW 87 0.11 0.05 0.84 Can-Seed 29 Can 17 EH 35 0.09 0.09 0.83 Can 15 Can-Sap-Seed 6 AB 98 0.01 0.95 0.04 Sap-Can-Seed 60 Sap-Seed-Can 20 IW 116 0.05 0.86 0.09 Sap-Seed-Can 33 Sap-Seed 31 Negative WP 30 0.3 0.37 0.33 Seed 7 Sap-Seed 5 parabola QA 31 0.06 0.29 0.65 Can 14 Sap 5 PC 16 0.31 0.5 0.19 Sap 7 Seed 5 AE 47 0.51 0.23 0.26 Seed 17 Can 6 Positive SM 141 0.46 0.11 0.43 Seed-Can-Sap 51 Can-Seed-Sap 43 parabola RO 37 0.41 0.14 0.46 Seed 10 Can-Seed 9 RM 119 0.8 0.08 0.13 Seed-Can-Sap 35 Seed 31 Right WA 72 0.78 0.18 0.04 Seed 19 Seed-Sap 15 skewed CH 118 0.55 0.25 0.2 Seed 31 Sap-Seed 17 Table 4.4. Summaries of stand-specific relative abundance rankings by size class for 18 species surveyed in > 10 % of our 141 stands. Scientific name and authorship can be looked up by species code in Table 4.1. Species are grouped by the shape of their average relative abundance size class distribution, and n is the number of sites that a species was present on, in any size class. Seedling, sapling, and canopy % indicates the percentage of stands in which that size class had the highest relative abundance (seedling 0 – 137 cm tall, saplings > 137 cm and < 10 cm DBH, and canopy/mature, abbreviated as canopy, > 10 cm DBH). 1st and 2nd ranking refer to two most common rankings of relative abundances by size class (n is number of stands). The nomenclature presents the size classes in decreasing rank order (i.e., the first size category listed had the highest relative abundance), with size classes which had zero stems omitted. For example, 21 out of 77 stands in which balsam fir (BF) was present had the relative abundance ranking “Can-Sap-Seed” (canopy is highest relative abundance); as a left-skewed species, it has greatest relative abundance in the canopy, on average. 104 4.5 Discussion 4.5.1 Overview Our results suggest individual stands contain few of the total tree species present in northern hardwood forests for seedling, sapling, and canopy/mature size classes and underscore current regeneration concerns. Despite sampling 25 tree species across the landscape, most had limited geographic distribution and low abundance; fewer than ten species were sampled in most stands or achieved at least 1% average relative abundance in any size class averaged across sites (landscape scale). At the stand level, all three size classes averaged fewer than three effective common species (CS), and although SR was similar between seedlings and canopy/mature stems, natural species filtering of sub-canopy and disease-impacted species means canopy/mature species richness is likely to decline as these seedling cohorts recruit to the canopy. Diversity patterns furthermore appear to be self-perpetuating, as stands with lower canopy/mature diversity (CS of stems > 10 cm DBH, CS10) tended to have lower seedling and sapling diversity (Table 4.3). These results have concerning implications for future stand resilience of managed northern hardwood forests in Michigan. Our diversity metrics are similar, if not slightly lower, than past studies of northern hardwood forests. However, contextualizing our results in relation to past literature is hampered by inconsistencies with scales of analysis and diversity metric calculations and limited studies analyzing diversity in managed northern hardwood forests at comparable scales. Average stand- level tree species richness of our stands, for all size classes, was lower than values reported for managed and old growth stands in Crow et al. (2002), though Crow covered a smaller geographic range than analyzed here and reported values for aggregates of three stands, which may inflate species diversity metrics, if stand histories, site qualities, or other factors differ. Our 105 results suggest CS was generally lower than shown in Danyagri et al. (2019), although that comparison relies on modeled estimates rather than direct comparison of range and mean. Sugar maple in our stands dominated the seedling class similarly to the managed northern hardwoods characterized in Neuendorff et al. (2007) and Crow et al. (2002) but is less dominant in the sapling size class than described in these studies, both of which were conducted in a limited geographic region in the western Upper Peninsula of Michigan which has been shown to support higher densities of sugar maple (Henry et al., 2021). Rather, our saplings patterns are characterized by roughly even representation of sugar maple, ironwood, and American beech (Figure 4.1), similar to findings in Angers et al. (2005), near Ottawa, Canada, and Elenitsky et al. (2020), in north-eastern Michigan, both of which are in regions of relatively lower lake-effect snowfall and therefore likely higher deer browsing pressure. Our results add to the pool of studies suggesting that existing tree regeneration diversity is unlikely to yield diverse and resilient canopies across a broad geographic range. We found mixed support for our first hypothesis that tree diversity would be highest for canopy/mature stems and on low quality sites. Seedling and canopy/mature stems tied for highest SR, but saplings had highest CS (Figure 4.2). However, although saplings are generally more equally represented by different species (Figure 4.1), four of the six species of highest relative abundance have low likelihood of attaining our canopy size class (> 10 cm DBH) and are even less likely to reach the upper canopy stratum of northern hardwood forests (i.e., 25 - 35 m tall) due to size-dependent pest and pathogens (i.e., American beech, ash species) or small maximum height (i.e., subcanopy species ironwood and balsam fir). Furthermore, although regeneration SR was comparable to the canopy, the actual pool of species with the potential to reach the canopy is smaller as seedlings must successfully transition through the sapling class to the canopy class, 106 and transition to sapling class is a bottleneck driven by several factors, including deer browsing (this study, Walters et al 2022) and deep shade/high BA (this study, Henry et al 2021, Walters et al. 2022). Overall, patterns across size classes suggest a future decline in canopy species richness. Patterns of diversity with site quality matched our hypotheses in the EUP, with CS generally declining with increasing site quality, but not in the other regions and only significantly for seedling and canopy/mature CS. Failure to detect patterns of diversity with site quality in other regions may be due to associations between high-quality sites and other factors. For example, deer use was positively associated with site quality in the WUP (Figure 1 in Henry et al., 2021). However, this does little to explain the lack of a negative site quality-sapling diversity relationship in the WUP as both high deer use and high site quality would be expected negatively associated with sapling diversity which we did not see. However, expectations that deer diminish diversity are not supported (this study). Via selective browsing deer likely impact composition more than diversity (this study; Walters et al., 2022). Our results support our overall second hypothesis, that tree species importance varies by size class (Figure 4.1), and many of our proposed mechanisms consistent with specific size pattern shapes were supported as well (Hypotheses 2a-d). Regarding Hypothesis 2a, we find that local seed source limitations may play a dominant role in driving regeneration composition and density. Longer-distance dispersal of seed from adjacent stands is unlikely, particularly for large- seeded, passively dispersed species such that local/intra-stand seed sources likely dominate seed availability (McEuan & Curran, 2002; Willis et al, 2016), and certainly sprout availability. Across all species and most regeneration size classes, we found positive associations with conspecific canopy/mature density for both large and small seeded species, and overall positive associations of seedling and sapling diversity measures with canopy/mature species diversity. 107 This finding supports the notion that local canopy seed/sprout sources have a strong bearing on seedling/sapling composition. Fewer significant associations between SR10 and regeneration diversity metrics as compared to CS10 suggests that canopy/mature density, rather than presence alone, is key for driving diverse regeneration. A caveat to our local seed/sprout source limitation interpretation is that other factors could contribute to canopy-sapling-seedling compositional continuity including site factors (e.g., soil water regimes) that favor establishment/growth/survival from seedling to maturity. (Webster et al., 2018). Our results, particularly our relative abundance analyses, suggest substrate limitation for small-seeded species, supporting our hypothesis 2b. Many small-seeded species, including eastern hemlock, white spruce, and birch species, had left-skewed relative abundances, with lower relative abundance in the understory compared to the overstory, supporting potential substrate limitation. Some of those species, such as eastern hemlock, have well-documented narrow seedling establishment substrate requirements with those substrates currently uncommon in managed forests (summarized in Alverson et al., 2019), which place our results in accordance with existing literature. However, distribution of relative abundance by size class is likely driven by multiple interacting factors and must be interpreted carefully. In addition to seedling establishment limitations, other species traits/properties may also contribute to left-skewed distributions. For example, for balsam fir and spruce high long-term survival, low growth rates for saplings and low deer browsing pressure may lead to the accumulation relatively high sapling densities. Another example, many of our small-seeded species (which is generally associated with narrow establishment substrate requirements: Marx and Walters, 2008; Bolton and D’Amato, 2011; 108 Willis et al., 2016) are also highly shade intolerant (e.g., PB, BTA) such that seedling and sapling densities are constrained by both establishment substrate and high low light mortality. Although substrate-limitation for small seeded species was supported by several relative abundance patterns, we identified few significant associates with species regeneration abundance and forest substrate composition (percent hardwood litter (HWL) or bare mineral soil (BMS)). However, several factors may explain this, particularly overall very low average coverage of BMS and CWD (Table 4.1). Less common species also had fewer non-zero data points, making association more difficult to detect, and our broad height range of seedlings (< 137 cm tall) may include stems that established decades ago (e.g., stems can persist < 1 m for > 30 years; Marks and Gardescu, 1998), which would potentially only weakly (or not at all) associate with current substrate conditions. Lastly, substrate limitation may be difficult to detect at a stand-scale, instead better analyzed at the plot level, as in Willis et al. (2016). Nevertheless, patterns of relative abundance suggest potential for seedling establishment to prove limiting in our stands dominated by hardwood litter coverage. Consistent with hypothesis 2c, high BA (i.e., low light) may strongly limit the transition from seedling to saplings for mid-shade tolerant species, including cherry, red maple, red oak, and tolerant sugar maple. Moderate to high shade tolerance can facilitate development of multiple, longer-lived seedling cohorts, which can accumulate as dense seedling populations; however, tolerance is not sufficient to support survival and growth into the sapling class, at least not in the low light environments characterizing dominant single tree selection silviculture in Great Lakes NHF. In general, negative associations of density with BA, positive associations with BASD and positive parabolic or right-skewed relative abundance distributions for mid- tolerant cherries, maples, and red oak all support the notion of low light limitation of seedling- 109 sapling recruitment. In addition to direct low-light/high BA limitations (i.e., mortality via carbon starvation) to seedling-sapling transitions for mid tolerant species, low-light conditions could exacerbate resource shortfalls for tree seedlings by negatively impacting arbuscular mycorrhizal (AM) colonization (Cheng et al., 2005; Koorem et al., 2017; Neuenkamp et al., 2021); red maple, cherry species, and sugar maple are all AM associated species (Bennett et al., 2017). In addition to high BA (light) limitations affecting relative abundance patterns across size classes, deer browsing may also drive declines in relative abundance from seedling to saplings for maples, cherries, and red oak, partially supporting hypothesis 2d. This pattern is most clearly supported for sugar maple, given significant negative associations between sapling density and deer use. Previous studies suggest deer are a likely contributing factor in regeneration bottlenecks, including for sugar maple (Côté et al., 2004; Leak, 2006; Neuendorff et al., 2007; Powers & Nagel, 2009; Matonis et al., 2011) oak species (Dey, 2014) and species groups including the maples, cherries and red oak examined here (Walters et al. 2022). In contrast to the positive parabola and right skewed forms for maples, cherries and red oak, several species know to be little impacted by deer browsing pressure (ironwood, American beech, white pine) had positive parabolic shape among size classes (i.e., saplings most abundant) which also supports the notion that deer are important drivers of compositional shifts. It is important to note that our species abundance model results regarding deer effects provided only limited support for our interpretation of deer effects based on size class changes in relative abundance. Sugar maple was the only species to negatively associated with DEER, and surprisingly, red oak saplings had a positive association. Our modeling approach may have failed to detect associations, in part, due to low stand-level detection and low abundances of many species (e.g., yellow birch). Densities 110 are also influenced by other factors previously discussed, such as low light and substrate limitation, which can further mask the effect of deer browsing. 4.5.2 Management implications In the context of forest management, our results support several potential alternative management strategies aimed at overcoming seed, substrate, and growth/survival (e.g., BA, deer browsing) limitations to establish species-diverse seedling and sapling classes. Given the current paucity of many tree species in all strata (and critically in the seed producing mature/canopy class), active seeding or planting of additional tree species may be necessary to promote greater diversity. Tree species selection for planting could be coupled with assisted migration strategies to further prepare forests for a variety of future climate scenarios. Provisions for exposing more mineral soil during harvest operations and leaving more trees to die and become coarse woody debris on the forest floor would help overcome establishment substrate limitations for several species, especially those with small seeds (Marx & Walters, 2008; Willis et al., 2015). Reducing basal area more than is done with current partial harvesting systems in stands would result in higher light levels than currently created following harvest and possibly increase sapling recruitment of shade mid-tolerant and tolerant species. Finally, implementing alternative strategies to deter deer browsing, such as leaving treetops as natural barriers, may increase regeneration density and recruitment of species already present on the landscape. Integrated strategies including several of the management modifications are likely to be most effective in increasing species diversity in NHF in the long term. Another management implication of our results is that the high variability we found across species in seedling vs sapling abundances (and ultimately to canopy/mature classes) emphasizes the potential risk in assuming seedling composition reflects sapling composition. 111 This is particularly important for developing meaningful assessments of regeneration success following harvest. Several factors, including shrub competition (Walters et al. 2016), deer browsing (Henry et al. 2021, this paper) and light resource limitation (i.e., canopy BA, Henry et al. 2021, Walters et al. 2022, this paper), may assert strong and variable bottlenecks on seedling to sapling transitions, such that tree seedlings are poor predictors of regeneration success (Walters et al. 2020a). Managers should focus on saplings (0 – 10 cm DBH) for assessment of regeneration success following harvest or other management activity. With certain key species potentially in decline, relative abundance analyses can be used generally to identify species occupying their relinquished growing space, namely right-skewed species (increasing relative abundance with decreasing size class) and negative parabolic species (with saplings highest relative abundance). Right-skewed species include red maple, which has been increasing in relative abundance in recent decades (Fei & Steiner, 2007). Many of the species which compete with sugar maple for sapling growing space are negative parabolic, including American beech, which creates well-documented dense thickets (Gravel et al., 2011; Collin et al., 2017; Elenitsky et al., 2020), and ironwood, which is a sub-canopy species. Relative abundance analyses prove to be a useful method of quantifying potential species trajectories, when incorporated with literature and knowledge of tree species’ life history dynamics. 4.5.3 Caveats Our snapshot-in-time survey is unable to assess trends within size classes over time, and instead relies on assumptions about current patterns and implications for ultimate recruitment; previous studies have suggested comparing size classes from a single survey may not be appropriate for inferring temporal compositional shifts (in the case of tree range shifts, Malis et al., 2016). With this caveat, our results nonetheless suggest that current managed northern 112 hardwood forests may face long-term challenges with resilience due to low diversity of species without proactive management interventions. Furthermore, our analysis is explicitly on managed systems, where seedling and sapling classes provide the raw material for management options aimed at future compositional goals, which enhances the utility/relevance of such analysis. However, follow-up analyses of species diversity dynamics over time would be useful in confirming the temporal component suggested by our relative abundance analyses, given the long-term maintenance of current compositional drivers over time (e.g., current management regimes, deer use). Nevertheless, our results highlight the importance of intervention to promote diverse, resilient stands which can adapt to changing timber market shifts, pest/pathogen introductions, and climatic changes. Applying ecological knowledge to these stands is vital for continued ecological functioning in managed northern hardwood forests. 113 CHAPTER 5 CONCLUSIONS The preceding chapters explored the outcomes of 60 + years of single-tree selection management, across a broad geographic gradient, in driving forest regeneration and recruitment outcomes. Chapters 2 and 3 emphasized that regeneration and recruitment have repeatedly failed for the focal species, sugar maple. Not only is sugar maple poorly regenerating, but it is being outcompeted by species less desirable for management which often do not or cannot recruit to the canopy. Overall, stands are characterized by low tree species diversity (Chapter 4). These results have critical implications for northern hardwood forest management. Analysis of sugar maple regeneration patterns and associations with key drivers confirms that sugar maple regeneration failure is geographically widespread and occurs on a variety of site conditions. While densities of sugar maple seedlings (< 50 cm tall) were generally high (though potentially limited on some sites by low mature sugar maple seed sources and/or competing taller vegetation), sugar maple stems > 50 cm tall and < 5 cm DBH were sparse and often entirely absent. Densities were particularly sparse in portions of the northern Lower Peninsula and the central southern Upper Peninsula; they also associated negatively with increasing winter deer use and understory competition, and positively with seed availability and medium site qualities. Our results echo previous findings of low sugar maple regeneration densities (Leak, 2006; Neuendorff et al., 2007; Powers & Nagel, 2009; Matonis et al., 2011), although differed from (Vickers et al., 2019) which analyzed all northern hardwood forests rather than being focused on selection managed stands. Our results emphasized that current stands are poorly stocked and unlikely to successfully regenerate new cohorts of sugar maple, given current stand dynamics. 114 Analysis of age structure of 1499 sugar maple basal discs from 51 stands indicated that single-tree selection has chronically failed to recruit new age cohorts of sugar maple in the past 60 + years of management (Chapter 3). As assumed from Michigan forest and logging history, most stand canopies are dominated by stems 90 – 120 years old; however, in contrast with expectations under uneven-aged single-tree selection management, we identified few age cohorts < 70 years old (36 of 160 identified age cohorts), with no stands stocked in stems 30 - 50 years old and only 7 of 51 stands stocked in stems 50 – 70 years old. Saplings average 66 years old, poletimber averages 91 years, and sawtimber averages 106 years. Regions with greater January precipitation (snowfall) and lower annual temperatures, where deer populations have been historically low, generally have younger saplings. These results contradict previous findings of sugar maple age structure in managed northern hardwood forests of the Great Lakes region, which typically identified young (< 60-year-old) saplings and linearly related age and diameter (Tubbs, 1977a; Dey et al., 2017; Harmala, 2021). Our results indicate that sugar maple recruitment has persistently failed in northern hardwood forests managed by single-tree selection, and relying on diameter distribution to assess regeneration and recruitment success may often be misleading; our results have concerning implications for future growth of aging northern hardwood stands. Despite declines in sugar maple understory dominance, few desirable species appear to have taken its place (Chapter 4). After accounting for sample size, stand-level species richness was slightly lower for saplings (> 137 m tall and < 10 cm DBH) versus seedlings (< 137 m tall), though both classes were comparable to canopy/mature species richness (> 10 cm DBH), and all size classes generally averaged fewer than three effective common species. Stands with lower canopy species diversity tended to have lower regeneration species diversity, suggesting 115 diversity is self-perpetuating. Our results suggest likely influence of substrate and light limitation, given patterns of relative abundance by size class for individual species and secondarily model results for individual species abundance. At the stand level, sugar maple demonstrates a “bottleneck” of lower sapling diversity compared to seedlings or canopy/mature stems, potentially driven by a legacy of deer browsing pressure and low-light levels prohibiting gap capture. These dynamics are similar to those found in past studies of managed northern hardwood forests (Angers et al., 2005; Danyagri et al., 2019; Elenitsky et al., 2020), and diversity is generally lower than past analyses of old growth northern hardwood forests (Crow et al., 2002). Should regeneration diversity patterns sustain through canopy recruitment, they signal future shifts in mature canopy stratum composition. Our results highlight a need for changes in management aimed at increasing tree species diversity. Together, these results suggest that continued selection silviculture in many areas would likely result in continued decline in species diversity and potential canopy failure, due to low density of desirable stems for management. Across our analyses, light limitation is a consistent associated factor with undesirable regeneration and recruitment outcomes. This suggests that stands may benefit from more intense harvesting which results in greater light availability. Based on our results, this could benefit not only increased recruitment of shade-intolerant species but the relatively shade-tolerant sugar maple species as well. Deer browsing also proves to be a consistent threat to northern hardwood forest regeneration. Since direct reduction of the deer population is socially unlikely, silvicultural prescriptions intended to reduce the impacts of deer on regeneration, such as treetops acting as natural browsing barriers, may improve regeneration outcomes. Lastly, the low diversity of canopy/mature tree species and associated lack of seed source appears to limit both large- and small-seeded species. Given lack of local seed source, 116 future management would benefit from direct seeding or planting to increase future forest diversity; this could also benefit future climate adaptations by selecting species which are projected to be better suited to the future climate, i.e., assisted migration. However, these efforts are likely to be challenged by continued, persistent deer browsing. Although there are numerous potential directions for related research, two avenues emerged as particularly promising areas for future work. First, given the unprecedented results of our sugar maple age structure analysis, additional research quantifying ring width of the samples would be prudent. This could elucidate changes in sugar maple growth and recruitment over time and could further explore past stand dynamics, such as growth response to disturbance (harvest). The second area of research is review and quantification of species diversity goals for northern hardwood forests. Although studies of northern hardwood forests and discussions of the importance of species diversity proliferate, there are no clear guidelines on what constitutes adequate species diversity to promote forest resilience. Furthermore, existing literature lacks consistency in past reporting of species diversity measures, which hinders efforts to define successful levels of tree diversity. These areas of research would further our understanding of forest stand dynamics in managed northern hardwood forests and guide future management. In conclusion, these results emphasize that single-tree selection management has not resulted in desired outcomes. Namely, stands have low tree diversity and low densities of desirable sugar maple stems, which have persistently failed to recruit to the canopy. Forest managers in many regions will likely need to consider alternative silvicultural regimes in the immediate future to sustain forest health and functioning. 117 APPENDICES 118 APPENDIX A Details on the deer model 119 A.1 Modeling framework We developed a two-stage modeling approach: Stage 1 was a deer model used to impute pre-harvest deer occurrence information to all stands; Stage 2, was a set of sugar maple regeneration models that were conditioned on Stage 1 (deer model) and a set of other key predictor variables. More specifically, given pre-harvest deer use surveys for a subset of our stands (a key sugar maple regeneration predictor), the Stage 1 stand-level deer use model provided point estimates of deer use, as predicted by climate and landcover variables, post- harvest deer use measurements, and a spatial random effect. Stage 2 comprised a multilevel model for plot-level sugar maple regeneration counts at three sizes classes as predicted by stand- level deer use from Stage 1 and stand vegetation structure and site quality measurements. All analyses and modeling were conducted in R (R Core Team, 2018). With a Bayesian framework, we modeled deer use prior to timber harvest (with deer use defined as the percent of transect segments with deer fecal pellets) as predicted by climate, landcover, post-harvest deer use, and a spatial random effect. Given the purpose of this model was to infer pre-harvest deer use in non-surveyed stands, we aimed to build the best predictive model rather than explain biological processes. Modeling percent of occupied transect segments (i.e., binary responses for each transect segment and rate or percent when considering multiple transect segments within a site) rather than counts or “deer use days,” the latter of which is common in the literature (e.g., (Yañez-Arenas et al., 2012), offers several advantages. Using simple pellet occurrence within a transect segment removes the need to count individual pellets or judge what constitutes a pellet group while incorporating uniformity of deposition. Moving to pellet occurrence also mitigated the potential for unusually abundant pellet group counts in a single segment which tended to provide unrealistically high stand-level deer density estimates. 120 This approach remains directly comparable to previous models based on counts (linear relationship between percent and counts for our data, R2 = 0.90). The bounded response variable expressed as a percent or rate also provides straightforward interpretation. We used a spatially stratified sampling approach (Fig. A.1). Despite a relatively broad range of sample dates, we conducted linear regression to test whether sampling date correlated with total pellet group counts and found no relationship (data not shown). Figure A.1. Design of the winter pellet survey transects. Points represent 25 sampling locations for vegetation surveys within the 12.14 ha square unit. Each transect (blue line) was 6 m wide, and cumulative transect length was 628 m. For analysis, transects were divided into spatially referenced segments ~ 6 m long. Occurrence of deer pellets were documented for each segment. The predictors for pre-harvest (2017) deer use include average daily snow depth (01 November 2016 - 30 April 2017, 1 km resolution, (National Operational Hydrologic Remote Sensing Center, 2004); average minimum monthly winter temperature (December 2016 - February 2017, 4 km resolution, 1981-2010, PRISM Climate Group, Oregon State University, http://prism.oregonstate.edu, created 18 October 2019); 2019 (post-harvest) deer use by timber harvest type; and local coverage of deciduous, evergreen, and mixed forests as well as agricultural lands within 25 km2 of the stand (NLCD, 30m resolution, aggregated to 1 km resolution for faster processing). These are well-documented drivers of deer populations in our study region (snow depth (Moen, 1976; Shi et al., 2006); winter low temperatures (Holter et al., 121 1975); and proximity to coniferous, deciduous, and agricultural land (Telfer, 1970; Ozoga & Gysel, 1972; Wetzel et al., 1975; Millington et al., 2010; Dawe et al., 2014). Although we acknowledge that local deer use can be concentrated by timber harvest activities (e.g., (Campbell et al., 2004), post-harvest deer use of a stand reflects relative availability of deer occupying the surrounding landscape thereby serving as a useful proxy for browse pressure. For the complete collection of predictor variables, we calculated variance inflation factors using the vif function in the car package (Fox & Weisberg, 2011) to check for multicollinearity; at a threshold of 2, we removed minimum winter temperature. With the remaining predictors, we modeled pre-harvest deer use (percentage of transect segments occupied by deer pellets at each stand, 𝑦𝑑𝑒𝑒𝑟,𝑖 ) using a polya-gamma binomial regression (Polson et al., 2013) and a spatially varying random intercept (Finley & Banerjee, 2020). The model for the ith stand was: 𝑦𝑑𝑒𝑒𝑟,𝑖 𝑃𝐺(𝑛𝑖 , 𝜓𝑖 ) 𝑙𝑜𝑔(𝜓𝑖 ) = 𝜇𝑑𝑒𝑒𝑟 + 𝛽𝑑𝑒𝑒𝑟,𝑖 𝑋𝑑𝑒𝑒𝑟,𝑖 + 𝑤𝑖 where 𝑛𝑖 is the number of trails/transect segments (𝑛𝑖 is 103 for all stands) and 𝜓𝑖 is the stand- specific likelihood of any transect segment being occupied (success); the log of 𝜓𝑖 is defined by a linear equation including the grand mean 𝜇𝑑𝑒𝑒𝑟 , vector of regression coefficients 𝛽𝑑𝑒𝑒𝑟,𝑖 and associated predictor design matrix 𝑋𝑑𝑒𝑒𝑟,𝑖 , and the spatially varying random intercept 𝑤𝑖 . Predictive inference was based on 50,000 post burn-in samples from three Markov chain Monte Carlo (MCMC) chains. We assessed model convergence using Gelman Rubin diagnostics in the coda package (Plummer et al., 2006), trace plots, and plots of residuals vs. predicted, all of which indicated model fit. Post burn-in samples were then used to generate predicted 𝜓𝑖 for all stands for subsequent use in the Stage 2 vegetation model. 122 A.2 Results and conclusions Across most of the study area, winter deer use during the winter of 2017 was estimated as low to moderate, except for relatively high use in the south-central and eastern parts of the Upper Peninsula (Fig. 2). Our model estimated an effective spatial range of 7.7 km (95% CI 2.8 – 27.5 km), indicating that sites > 7.7 km apart on average were not spatially autocorrelated (Table A.1). None of the weather or landcover variables significantly predicted deer use, and only 2019 deer pellets in single-tree selection harvest were positively associated with 2017 pellet counts (Table A.1). Our deer use predictions (Fig. 2) aligned with other smaller-scale studies of deer use in our study region (Shi et al., 2006; Millington et al., 2010); this, along with reasonable estimates of spatial autocorrelation (~ 7 km) that correspond to deer winter habitat range (Van Deelen et al., 1998), gave us confidence that our model reasonably predicted deer use for sites in our study area. 123 Predictor Median (95% CI) Intercept -1.013 (-2.590 – 0.663) Average snow depth -0.015 (-0.038 - 0.008) Deciduous forest cover -0.915 (-2.606– 0.739) Evergreen forest cover -1.129 (-6.265 – 3.594) Mixed forest cover 4.407 (-1.354 – 10.520) Agricultural land cover 3.492 (-4.621 – 12.052) 2019 pellets (single-tree selection harvest) 0.019 (0.000 – 0.037) 2019 pellets (seed tree harvest) 0.007 (-0.007 – 0.020) 2019 pellets (shelterwood harvest) 0.011 (-0.003 – 0.024) 2019 pellets (group selection harvest) 0.013 (-0.002 – 0.029) Sigma squared 0.695 (0.441 – 1.146) Effective special range (km) 7.086 (3.789 – 24.984) Table A.1. Results from the stage one deer use model are presented with their median value and a 95% Bayesian confidence interval. Sigma squared is the variance parameter for the spatial term, and effective spatial range represents the distance at which spatial autocorrelation becomes insignificant in the model, in kilometers. None of the climatic or landscape predictor variables were significantly different from zero, likely because the spatial term explained most of the variation. 124 APPENDIX B Additional information on site quality categorizations 125 Site Quality Abbrev. NLP EUP WUP TOTAL Poor - Q1 PArVVb 4 PArVAa, 27 pArVAa, 8 39 poor/medium ATFD, AArAst, AArLy, PArVAa[w] Medium* Q2 - - AFPo 10 ATM-Sm, ATM- 12 22 O, ATFAs Medium/rich Q3 AFO 23 AFOAs 8 ATD-Ca, ATD- 16 47 Hp, AVb Rich - very rich Q4 AFOCa 21 - - AOCa 12 33 TOTAL 48 45 48 141 Table B.1. Bayesian multilevel model predicting sugar maple regeneration in Michigan included four ordinal categories of site quality (sensu (Burger & Kotar, 2003). Kotar site quality grouping, model abbreviation, Kotar types (see Table B.2) included in the grouping, and count of stands sampled by region. The full set of site qualities was only present in one of the regions. NLP = Northern Lower Peninsula, EUP = Eastern Upper Peninsula, WUP = Western Upper Peninsula. *Medium site quality was included in the intercept. 126 Habitat type Region Primary Primary soils Soil Soil landforms moisture nutrient regime regime PArVVb (Pinus strobus – NLP Outwash Moderately Dry- Poor to Acer rubrum/Vaccinium- plains, ground well drained Mesic medium Viburnum acerifolium) moraines, and sandy soils beach ridges AFO (Acer saccharum – NLP End moraines, Well drained Mesic Medium Fagus ground sandy soils to rich grandifolia/Osmorhiza moraines, and claytoni) outwash plains AFOCa (Acer saccharum NLP End and Well to Mesic Rich to – Fagus ground moderately well very rich grandifolia/Osmorhiza moraines drained deep claytoni-Caulophyllum sandy loam till thalictroides) PArVVa (Pinus strobus – EUP Deep Excessively to Dry to Poor Acer rubrum/Vaccinium lacustrine well drained dry- angustifolium – Aralia deposits of sandy soils mesic nudicaulis) sand and gravel ATFD (Acer saccharum – EUP Outwash, Well to Mesic Poor to Tsuga canadensis – lacustrine moderately well medium Fagus deposits, drained deep grandifolia/Dryopteris glacial till, sands spinulosa) and end morains AFPO (Acer saccharum – EUP Variety of Well to Mesic Medium Fagus landforms somewhat grandifolia/Polygonatum excessively pubescens) drained deep sands and loamy sands; gravelly, cemented, and mottled layers common Table B.2. Descriptions of site quality types (Burger and Kotar 2003) used to model sugar maple regeneration. NLP = northern Lower Peninsula, WUP = western Upper Peninsula, and EUP = eastern Upper Peninsula. 127 Table B.2 [cont’d] Habitat type Region Primary Primary soils Soil Soil landforms moisture nutrient regime regime AFOAs (Acer saccharum – EUP End moraines Moderately Mesic Medium Fagus and till plains well to to rich grandifolia/Osmorhiza excessively claytoni – Arisaema drained soils; atrorubens) gravelly, cemented, and mottled layers common; thin till over bedrock pArVAa (Pinus strobus – WUP Glaciofluvial Excessively Dry to Poor Acer rubrum/Vaccinium deposits, well drained dry- angustifolium – Aralia moraines, soils of mesic nudicaulis) lake plains lacustrine deposits AArAst (Acer saccharum – WUP Coarse till Sandy soils Dry- Poor to Acer rubrum/Aster and shallow mesic medium macrophyllus) till over bedrock AArLy (Acer saccharum – WUP Coarse till Loamy soils Mesic Poor to Acer rubrum/Lycopodium deposits and medium annotinum) thin till over bedrock PArVAa[w] (Pinus strobus WUP Glacial Sand to sandy Dry- Poor to – Acer rubrum/Vaccinium outwash and loam mesic medium angustifolium – Aralia moraines nudicaulis [Wisconsin variant]) ATM-Sm (Acer saccharum WUP End moraines Loamy sand Mesic Medium – Tsuga and outwash and sandy canadensis/Maianthemum covered loam soils canadense – Smilacina moraines racemose variant) ATM-O (Acer saccharum – WUP Clay and Sandy loam Mesic Medium Tsuga lacustrine soils canadensis/Maianthemum deposits canadense – Osmorhiza claytoni variant) 128 Table B.2 [cont’d] Habitat type Region Primary Primary soils Soil Soil landforms moisture nutrient regime regime ATFAs (Acer saccharum WUP Lacustrine Sandy soils Mesic Medium – Tusga canadensis – deposits of with subsurface Fagus sand and clayey, grandifolia/Arisaema glaciofluvial gravelly, or atrorubens) deposits cemented layers ATD-Ca (Acer WUP Clay deposits Loamy cap soils Mesic to Rich saccharum – Tsuga wet- canadensis/Dryopteris mesic spinulosa – Caulophyllum thalictroides variant) ATD-Hp (Acer WUP Medium Sandy soils Mesic Medium saccharum – Tsuga textured with subsurface to rich canadensis/Dryopteris glacial till clayey, spinulosa – Hepatica gravelly, or variant) cemented layers AVb (Acer WUP Medium Sandy loams Dry- Medium saccharum/Viburnum textured end mesic to rich acerifolium) moraines AOCa (Acer WUP Moraines and Well drained Mesic Rich to saccharum/Osmorhiza loess deposits loamy till and very rich claytoni – Caulophyllum loess thalictroides) 129 APPENDIX C Sugar maple model: parameter and model selection 130 C.1 Parameter selection We compared linear models of log transformed stand-level average sugar maple seedling counts (+ 0.001 to account for zeros) predicted by different measures of canopy sugar maple abundance. For predictors, we compared average basal area and stem count of sugar maple trees within the 12.6 m2 plots (25 plots per stand) greater than thresholds of 25, 30, 35, and 50 cm DBH. Basal area of sugar maple canopy trees greater than 25cm DBH had the greatest predictive power and was included in the final model (Table C.1). Diameter threshold (cm) Stem count adjusted R2 Basal area adjusted R2 25 0.08871 0.1171 30 0.09429 0.1128 35 0.1031 0.1068 40 0.06418 0.07021 Table C.1. Adjusted R2 values for a linear models of stand-level log average sugar maple seedling counts predicted by different measures of canopy sugar maple density. Predictors included stand-level average stem count and basal area of sugar maple trees within 6 m2 radius plots, at varying size thresholds. The basal area of sugar maple trees > 25 cm had the highest adjusted R2 model and was used as the predictor of seed availability in model used to predict sugar maple regeneration. Based on (Elenitsky et al., 2020), we assessed whether densities of sugar maple regeneration in three size classes were comparable within ordinal site quality categorization, across three geographic regions for which our habitat classification system is parameterized (Burger & Kotar, 2003). We conducted Wilcoxon paired t-tests with a Bonferroni correction (n=24) (Table C.2). Only two of the 24 t-test comparisons significantly differed from zero, meaning that for those two stands, densities were not comparable within the site qualities across the two regions being compared. This suggests that overall regeneration densities are comparable across regions in the same site quality categorization. We also conducted recursive partitioning to indicate the most parsimonious groupings, in this case of site-level sugar maple densities as 131 predicted by region and ordinal site quality (Figs. C.1). If clear trends exist that regeneration is first grouped by region (higher at the top of the tree) and then by site quality, this would indicate that stands are not comparable within the same site quality categorization across regions; we failed to detect clear evidence of this pattern. Collectively, these results indicated that ordinal site quality by region was an appropriate predictor variable in our sugar maple regeneration model. Size Class Habitat class NLP- EUP- WUP- EUP WUP NLP Seedlings Poor to poor/medium 1 0.0078 0.3879 Medium - 1 - Medium/rich 1 1 1 Rich to very rich - - 1 Small Poor to poor/medium 1 .03 1 saplings Medium - 1 - Medium/rich 1 1 1 Rich to very rich - - 1 Small Poor to poor/medium 1 1 1 saplings Medium - 1 - Medium/rich 1 1 1 Rich to very rich - - 1 Table C.2. Results of the pairwise Wilcoxon tests using a Bonferroni correction (n=24). Values indicate probability of the distributions of the pair being significantly different. Values significant at adjusted p-values of 0.05 are highlighted. 132 A. Figure C.1. Decision tree generated by ANOVA model recursive partitioning for total stand-level counts of sugar maple seedlings (C.1A), small saplings (C.1B), and large saplings (C.1C). Predictive inputs were unique identifiers for site quality (1st, 2nd, 3rd, 4th, with 1st being the highest quality and 4th being the lowest quality) by region (eastern Upper Peninsula - EUP, western Upper Peninsula - WUP, and northern Lower Peninsula - NLP); for example, 1st EUP indicates highest quality sites in the eastern Upper Peninsula. In the bottom boxes, count values indicate the predicted value (total count of sugar maple stems surveyed per stand) while percentage values indicate the percentage of observations (out of 141) in the node. 133 Figure C.1 [cont’d] B. C. 134 We used AIC comparison (Akaike, 1974) to assess if interaction terms improved a simple model of regeneration predicted by deer use and site quality (glm.nb from MASS package, (Venables & Ripley, 2002). AIC comparison indicated that an interaction term would benefit the large saplings model (∆AIC -2, no significant interaction terms), but not for small saplings or seedlings, so we chose not to include the interaction term for any model (data not shown). C.2 Model selection We compared five candidate models for each sugar maple response size class (Table C.3). The models included a null model, a model with only linear terms, a model with a quadratic terms (^2) for each continuous predictor, a model with quadratic terms on all continuous predictors which may limit light (all continuous other than deer use and shrub coverage), and a model with quadratic terms based on prior AIC exploration of the data. For the model informed by AIC exploration, we compared the fit of linear vs. quadratic models (which include a linear and second-order polynomial component) for each variable individually with the response, at the stand-level, using glm.nb from the MASS package; we included quadratic terms in the candidate model based on R2 comparison. For each of the candidate models for each size class, we calculated the DIC value based on 10,000 post burn-in samples from three MCMC chains using the RJAGS package (Plummer, 2018). 135 Seedlings Candidate Model Model hypothesis DIC 1 SM1 ~ a; Null model 22228 a~1 2 SM1 ~ a + tot2 + tot3 + sap + wc; Linear effects only 22063 a ~ SM10BA + standBA + deer + Nut1 + Nut3 + Nut4 3 SM1 ~ a + tot2 + tot2^2 + tot3 + tot3^2 + sap + Full quadratic effects 22039 sap^2 + wc + wc^2; a ~ SM10BA + SM10BA^2 + standBA + standBA^2 + deer + deer^2 + Nut1 + Nut3 + Nut4 4 SM1 ~ a + tot2 + tot2^2 + tot3 + tot3^2 + sap + Quadratic terms for light 22030 sap^2 + wc; limitation a ~ SM10BA + SM10BA^2 + standBA + standBA^2 + deer + Nut1 + Nut3 + Nut4 5 SM1 ~ a + tot2 + tot2^2 + tot3 + tot3^2 + sap + Quadratic terms based on 22036 sap^2 + wc + wc^2; prior AIC exploration a ~ SM10BA + SM10BA^2 + standBA + deer + deer^2 + Nut1 + Nut3 + Nut4 Table C.3. Candidate models for sugar maple seedlings (SM1), small saplings (SM2), and large saplings (SM3) predicted by plot level densities of total small sapling stems (tot2), total large sapling stems (tot3), stems 5-10 cm DBH (sap), and shrub coverage (wc), plus the stand-level measurements of basal area of sugar maple stems > 25 cm DBH (SM10BA), total stand basal area > 10 cm (standBA), deer use (deer), and site quality (Nut1, Nut3, Nut4) which generate a stand-level intercept, a. The DIC values are based on 10,000 post burn-in samples from three MCMC chains using the RJAGS package (Plummer, 2018). For each size class, the model with the lowest DIC value is bolded. 136 Table C.3 [cont’d] Small saplings Candidate Model Model hypothesis DIC 1 SM2 ~ a; Null model 10620 a~1 2 SM2 ~ a + tot3 + sap + wc; Linear effects only 10590 a ~ SM10BA + standBA + deer + Nut1 + Nut3 + Nut4 3 SM2 ~ a + tot3 + tot3^2 + sap + sap^2 + wc + Full quadratic effects 10436 wc^2; a ~ SM10BA + SM10BA^2 + standBA + standBA^2 + deer + deer^2 + Nut1 + Nut3 + Nut4 4 SM2 ~ a + tot3 + tot3^2 + sap + sap^2 + wc; Quadratic terms for light 10175 limitation a ~ SM10BA + SM10BA^2 + standBA + standBA^2 + deer + Nut1 + Nut3 + Nut4 5 SM1 ~ a + tot3 + sap + wc; Quadratic terms based on 10486 prior AIC exploration a ~ SM10BA + SM10BA^2 + standBA + deer + Nut1 + Nut3 + Nut4 Large Saplings Candidate Model Model hypothesis DIC 1 SM3 ~ a; Null model 8676 a~1 2 SM3 ~ a + sap + wc; Linear effects only 8601 a ~ SM10BA + standBA + deer + Nut1 + Nut3 + Nut4 3 SM3 ~ a + sap + sap^2 + wc + wc^2; Full quadratic effects 8589 a ~ SM10BA + SM10BA^2 + standBA + standBA^2 + deer + deer^2 + Nut1 + Nut3 + Nut4 4 SM3 ~ a + sap + sap^2 + wc; Quadratic terms only on 8575 a ~ SM10BA + SM10BA^2 + standBA + light limiting variables standBA^2 + deer + Nut1 + Nut3 + Nut4 5 SM1 ~ a + sap + wc + wc^2; Quadratic terms based on 8611 a ~ SM10BA + standBA + deer + Nut1 + Nut3 + prior AIC exploration Nut4 137 APPENDIX D Details on species diversity analysis 138 Site quality WUP EUP NLP 1 AArAst, AArLy, PArVAa[w] (n=3) ATFD (n=10) PArVVb (n=2) 2 ATM, ATM-Sm, ATM-O, ATFAs (n=4) AFPo (n=5) - 3 ATD, ATD-Hp, ADT-Ca, ATD-AVb (n=6) AFOAs (n=3) AFO (n=8) 4 AOCa (n=4) - AFOCa (n=6) Table D.1. Summary of stand (n = 51) classification under an herbaceous-indicator site-quality system (Burger & Kotar, 2003). Site quality 1 is lowest, 4 is highest. Letters represent key indicator species: A = Acer saccharum, Aa = Aralia nudicaulis, Ar = Acer rubrum, As = Arisaema atrorubens, Ast = Aster marophyllus, Ca = Caulophyllum thalictroides, D = Dryopteris spinulosa, F = Fagus grandifolia, Hp = Hepatica [variant], Ly = Lycopodium annotinum, M = Maianthemum canadense, O = Osmorhiza claytoni [variant], P = Pinus strobus, Po = Polygonatum pubescens, Sm = Smilacina racemose [variant], T = Tsuga canadensis, V = Vaccinium angustifolium, Vb = Viburnum acerifolium, and [w] = Wisconsin variant. 139 Plot Scientific name Common name occurrences Abies balsamea Balsam fir 13 Acer rubrum Red maple 25 Acer saccharum Sugar maple 51 Betula alleghaniensis Yellow birch 13 Betula papyrifera Paper birch 2 Fagus grandifolia American beech 26 Fraxinus americana White ash 5 Fraxinus pennsylvanica Green ash 2 Juniperus virginiana Eastern red cedar 1 Ostrya virginiana Ironwood 10 Picea glauca White spruce 2 Pinus strobus White pine 1 Populus grandidentata Bit-tooth aspen 2 Populus tremuloides Quaking aspen 2 Prunus serotina Black cherry 10 Prunus spp. Cherry species 2 Quercus rubra Northern red oak 4 Tilia americana Basswood 25 Tsuga canadensis Eastern hemlock 2 Ulmus americana American elm 1 Table D.2. Summary of tree species (scientific and common name plus number of plot occurrences (out of 51)) from plot-level data, surveyed in 51 managed northern hardwood forests. 140 Stand Linear AIC CR AIC Pseudo R2 Stand Linear AIC CR AIC Pseudo R2 EUP1 196.9 187 0.806 NLP27 320.4 322.7 0.267 EUP2 228.5 226.7 0.417 NLP28 242.3 205.2 0.91 EUP3 236.7 239.8 0.332 NLP29 261.3 259.7 0.541 EUP4 210.4 207.7 0.429 NLP30 157.9 156.3 0.461 EUP5 233.7 232.7 0.482 NLP31 223.7 202.4 0.838 EUP6 267.3 265.8 0.49 NLP32 218.4 215.1 0.384 EUP7 218.6 204 0.645 NLP33 129.4 131.3 0.00101 EUP8 195.9 193.5 0.808 NLP34 243.8 232.6 0.687 EUP9 235.6 233.4 0.647 WUP35 345.5 348.1 0.345 EUP10 219.8 205.4 0.75 WUP36 236.3 235.8 0.631 EUP11 253.3 246.8 0.662 WUP37 327.5 328.5 0.658 EUP12 251.4 255 0.691 WUP38 304.5 306.9 0.624 EUP13 200.9 191.7 0.789 WUP39 247.2 234.9 0.744 EUP14 254 256.6 0.426 WUP40 245.5 235.6 0.642 EUP15 196.9 195.6 0.404 WUP41 195.1 167.7 0.861 EUP16 366.4 369.4 0.512 WUP42 267.6 266.7 0.427 EUP17 242.9 240.5 0.594 WUP43 225.7 209.2 0.698 EUP18 278 278.8 0.824 WUP44 176 154.4 0.881 NLP19 177.8 177.6 0.512 WUP45 271.4 272.8 0.454 NLP20 212.6 213.9 0.4 WUP46 289.2 288.8 0.224 NLP21 215.7 210.8 0.494 WUP47 251.8 222.8 0.785 NLP22 111.2 112.6 0.166 WUP48 246.5 210.8 0.922 NLP23 195.5 181.5 0.824 WUP49 250.3 245 0.671 NLP24 251.6 242.2 0.636 WUP50 229.5 229.7 0.182 NLP25 238.5 233.9 0.312 WUP51 303.4 307.6 0.738 NLP26 233.4 225.7 0.741 Table D.3. AIC comparison between linear and Chapman-Richards (CR) models, by stand, plus the pseudo-R2 value of the selected model (lower AIC, selected model AIC is bolded). Stands are identified by region (EUP = eastern Upper Peninsula, WUP = western Upper Peninsula, and NLP = northern Lower Peninsula). 141 Asym b c Stand E SE p-val E SE p-val E SE p-val EUP1 105.97 5.73 <0.001 0.1 0.05 0.042 0.96 0.49 0.061 EUP2 104.14 7.33 <0.001 0.18 0.13 0.182 1.66 1.93 0.400 EUP4 95.77 7.8 <0.001 0.08 0.08 0.357 0.43 0.39 0.273 EUP5 107.11 10.1 <0.001 0.09 0.07 0.205 1.12 1.09 0.310 EUP6 112.42 16.92 <0.001 0.07 0.07 0.330 0.67 0.51 0.200 EUP7 93.3 2.36 <0.001 0.24 0.09 0.016 1.58 1.12 0.171 EUP8 157.08 56.61 0.012 0.03 0.03 0.376 0.77 0.31 0.023 EUP9 138.74 26.38 <0.001 0.05 0.05 0.312 0.72 0.41 0.094 EUP10 105.57 4.19 <0.001 0.16 0.05 0.002 3.2 1.82 0.092 EUP11 128.29 7.36 <0.001 0.1 0.05 0.042 1.13 0.63 0.086 EUP13 106.76 4.61 <0.001 0.09 0.04 0.033 0.6 0.25 0.024 EUP15 125.77 4.2 <0.001 0.08 0.05 0.145 0.53 0.47 0.271 EUP17 140.31 18.31 <0.001 0.05 0.05 0.304 0.55 0.36 0.134 NLP19 128.67 20.56 <0.001 0.07 0.07 0.317 1.02 1.02 0.332 NLP21 97.84 3.79 <0.001 0.14 0.1 0.144 0.64 0.55 0.257 NLP23 110.73 3.77 <0.001 0.09 0.04 0.044 0.46 0.18 0.018 NLP24 98.09 3.65 <0.001 0.11 0.05 0.057 0.87 0.55 0.125 NLP25 102.25 3.08 <0.001 0.24 0.2 0.246 1.63 3.08 0.601 NLP26 116.38 4.79 <0.001 0.13 0.06 0.028 1.08 0.57 0.068 NLP28 103.44 2.29 <0.001 0.14 0.02 <0.001 2.09 0.55 <0.001 NLP29 102.31 12.97 <0.001 0.06 0.07 0.387 0.46 0.32 0.164 NLP30 106.31 6.38 <0.001 0.06 0.07 0.387 0.3 0.28 0.301 NLP31 100.03 2.99 <0.001 0.08 0.03 0.004 0.67 0.22 0.005 NLP32 101.94 2.9 <0.001 0.25 0.14 0.097 1.81 2.44 0.464 NLP34 108.04 4.43 <0.001 0.14 0.05 0.007 1.34 0.69 0.061 Table D.4. Model parameter estimates, including mean (E), standard error (SE), and p-value (p- val), for AIC-selected Chapman Richards models. Model estimates significantly different from zero (p < 0.05) are bolded. The model equation for age of stem i as a function of DBH is: 𝐴𝑔𝑒𝑖 ~ 𝐴𝑠𝑦𝑚 ∗ (1 − 𝑒 −𝑏∗𝐷𝐵𝐻𝑖 )𝑐 , where the variable Asym is interpretable as the age at which model estimates plateau. Stands are identified by region (EUP = eastern Upper Peninsula, WUP = western Upper Peninsula, and NLP = northern Lower Peninsula). 142 Table D.4 [cont’d] Asym b c Stand E SE p-val E SE p-val E SE p-val WUP36 96.1 26.12 <0.001 0.04 0.06 0.563 0.42 0.28 0.153 WUP39 96.05 4.44 <0.001 0.1 0.04 0.027 0.81 0.37 0.036 WUP40 96.67 2.85 <0.001 0.15 0.07 0.036 1.07 0.75 0.162 WUP41 99.22 1.1 <0.001 0.32 0.07 <0.001 4.41 2.74 0.121 WUP42 118.25 20.4 <0.001 0.1 0.08 0.220 2.33 2.52 0.365 WUP43 117.38 1.57 <0.001 0.4 0.16 0.017 15.54 24.36 0.529 WUP44 118.19 2.74 <0.001 0.12 0.03 0.001 1.34 0.48 0.011 WUP46 131.65 20.31 <0.001 0.09 0.11 0.428 0.9 1.2 0.459 WUP47 107.25 1.45 <0.001 0.22 0.06 <0.001 1.35 0.63 0.040 WUP48 107.72 3.44 <0.001 0.13 0.03 <0.001 1.93 0.49 <0.001 WUP49 103.71 6.82 <0.001 0.09 0.06 0.146 0.6 0.35 0.097 143 Intercept Slope Stand E SE p-val E SE p-val EUP3 36.3 7.97 <0.001 1.44 0.4 0.001 EUP12 35.83 7.01 <0.001 1.92 0.24 <0.001 EUP14 33.29 12 0.01 2.42 0.56 <0.001 EUP16 24.49 13.76 0.084 2.91 0.48 <0.001 EUP18 34.39 6.92 <0.001 2.97 0.25 <0.001 NLP20 80.47 5.17 <0.001 0.58 0.14 <0.001 NLP22 117.6 6.26 <0.001 0.31 0.17 0.083 NLP27 66.94 11.63 <0.001 1.53 0.45 0.002 NLP33 97.99 4.36 <0.001 0.02 0.14 0.891 WUP35 52.03 8.53 <0.001 1.63 0.37 <0.001 WUP37 36.74 12.82 0.007 4.3 0.56 <0.001 WUP38 30.5 10.77 0.008 2.58 0.37 <0.001 WUP45 56.18 7.8 <0.001 1.41 0.29 <0.001 WUP50 95.88 3.63 <0.001 0.4 0.15 0.015 WUP51 31.2 9.64 0.003 3.01 0.33 <0.001 Table D.5. Model parameter estimates, including mean (E), standard error (SE), and p-value (p- val) for AIC-selected linear models. Model estimates significantly different from zero (p < 0.05) are bolded. The model (which is unique for each stand) for age of stem i a function of DBH is: 𝐴𝑔𝑒𝑖 ~ 𝐼𝑛𝑡𝑒𝑟𝑐𝑒𝑝𝑡 + 𝑆𝑙𝑜𝑝𝑒 ∗ (𝐷𝐵𝐻𝑖 ). Stands are uniquely numbered and identified by region (EUP = eastern Upper Peninsula, WUP = western Upper Peninsula, and NLP = northern Lower Peninsula). 144 Variable Estimate Std. Error T value P value Intercept 80.2 2.5 31.5 < 2e-16 CPIntermediate 14.0 2.8 5.0 6.7e-07 CPCodominant 25.0 2.4 10.6 < 2e-16 CPDominant 42.7 2.9 14.7 < 2e-16 Table D.6. Model summaries of fixed effects of crown class as a factor predicting sugar maple age (overtopped included in the intercept, intermediate CPIntermediate, codominant CPCodominant, and dominant CPDominant stems). Significant predictor estimates (p < 0.05) are bolded. Variable Estimate Std. Error T value P value Intercept 74.64 6.00 12.5 < 2e-16 Crown percentage -0.31 0.13 -2.4 0.0183 DBHcm 0.97 0.21 4.6 5.48e-06 Interaction term 0.01 0.005 2.4 0.0190 Table D.7. Model summary of the fixed effect variables crown percentage, DBH (cm), and their interaction as predictors of sugar maple age. Significant predictor estimates (p < 0.05) are bolded. Variable Estimate Std. Error T value P value Intercept 55.7 3.79 14.73 < 2e-16 Height (m) 2.2 0.23 9.54 < 2e-16 Table D.8. Model summary of a fixed effects of sugar maple stem age as a function of height. Significant predictor estimates (p < 0.05) are bolded. 145 Figure D.1. Individual tree disc ages (n = 1499) by competitive class (O = overtopped, I = intermediate, C = codominant, and D=dominant). Figure D.2. Map of model-averaged sugar maple stem age at 10 cm DBH for 51 managed northern hardwood forests. 146 Figure D.3. Maps of average sample age for sugar maple stems in 51 managed northern hardwood forests in three size classes: saplings (5 – 11.4 cm DBH), poletimber (11.4 – 23.9 cm DBH), and sawtimber (> 23.9 cm DBH). 147 APPENDIX E Model selection for stand-level tree species abundance 148 DIC Species Size class NB ZINB AB 1 1193.660342 2225.652189 AB 2 1598.378344 2109.714133 AE 1 4.44771E+12 54769793063 AE 2 4.34462E+11 4027.4578 BF 1 17662.55378 5859.864774 BF 2 6956.414702 1540.43038 BW 1 1164.105206 1335.87357 BW 2 8948264249 8203.240323 CH 1 1379.812449 4515.791002 CH 2 1004.47792 2404.543577 IW 1 1305.503868 3883.088101 IW 2 1638.948521 3830.793844 RM 1 1829.801164 7294.77764 RM 2 1152.69937 4105.586483 RO 1 7.17E+25 562.488586 RO 2 2.48E+48 1.12107E+12 SM 1 2457.502827 6065.31646 SM 2 1912.313374 2517.543047 WA 1 43203411.13 3863.874594 WA 2 1.11E+20 6661.727899 YB 1 2455971.93 9078.715043 YB 2 162518909.4 6156.992703 Table E.1. DIC comparison for negative binomial (NB) and zero-inflated negative binomial (ZINB) models. Lower DIC value is bolded. Species codes can be found in Table 4.1. Size class 1 refers to seedlings (0 – 137 cm tall) and size class 2 refers to saplings (> 137 cm tall and < 10 cm DBH). 149 WORKS CITED 150 WORKS CITED Aitken, S. N., Yeaman, S., Holliday, J. A., Wang, T., & Curtis-McLane, S. (2008). Adaptation, migration or extirpation: climate change outcomes for tree populations. Evolutionary Applications, 1, 95–111. https://doi.org/10.1111/j.1752-4571.2007.00013.x Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6), 716–723. https://doi.org/10.1109/TAC.1974.1100705 Allen, C. D., Macalady, A. K., Chenchouni, H., Bachelet, D., McDowell, N., Vennetier, M., Kitzberger, T., Rigling, A., Breshears, D. D., Hogg, E. H. (Ted), Gonzalez, P., Fensham, R., Zhang, Z., Castro, J., Demidova, N., Lim, J. H., Allard, G., Running, S. W., Semerci, A., & Cobb, N. (2010). A global overview of drought and heat-induced tree mortality reveals emerging climate change risks for forests. Forest Ecology and Management, 259(4), 660–684. https://doi.org/10.1016/j.foreco.2009.09.001 Alverson, W. S., Lea, M. v., & Waller, D. M. (2019). A 20-year experiment on the effects of deer and hare on eastern hemlock regeneration. Canadian Journal of Forest Research, 49(11), 1329–1338. https://doi.org/10.1139/cjfr-2019-0071 Angers, V. A., Messier, C., Beaudet, M., & Leduc, A. (2005). Comparing composition and structure in old-growth and harvested (selection and diameter-limit cuts) northern hardwood stands in Quebec. Forest Ecology and Management, 217(2–3), 275–293. https://doi.org/10.1016/j.foreco.2005.06.008 Arbogast, C. (1957). Marking guides for northern hardwoods under the selection system. In USDA Forest Service Station Paper No. 56. Arseneault, J. E., Saunders, M. R., Seymour, R. S., & Wagner, R. G. (2011). First decadal response to treatment in a disturbance-based silviculture experiment in Maine. Forest Ecology and Management, 262(3), 404–412. https://doi.org/10.1016/j.foreco.2011.04.006 Bannon, K., Delagrange, S., Belanger, N., & Messier, C. (2015). American beech and sugar maple sapling relative abundance and growth are not modified by light availability following partial and total canopy disturbances. Canadian Journal of Forest Research, 45, 632–638. Baral, S. K., Danyagri, G., Girouard, M., Hébert, F., & Pelletier, G. (2016). Effects of suppression history on growth response and stem quality of extant northern hardwoods following partial harvests. Forest Ecology and Management, 372, 236–246. https://doi.org/10.1016/j.foreco.2016.04.023 Baribault, T. W., & Kobe, R. K. (2011). Neighbour interactions strengthen with increased soil resources in a northern hardwood forest. Journal of Ecology, 99(6), 1358–1372. 151 Baribault, T. W., Kobe, R. K., & Rothstein, D. E. (2010). Soil calcium, nitrogen, and water are correlated with aboveground net primary production in northern hardwood forests. Forest Ecology and Management, 260, 723–733. https://doi.org/10.1016/j.foreco.2010.05.029 Bartholomew, C. (2018). Overtopped, volunteer, 30- to 40-yr-old sugar maple saplings responded to unplanned release events in an even-aged, maturing sugar maple/northern red oak plantation. https://digitalcommons.esf.edu/etds Barton, K. (2020). MuMIn: Multi-Model Inference. R package version 1.43.17. https://CRAN.R- project.org/package=MuMIn. Bates, D., Maechler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48. Beaudet, M., Messier, C., & Canham, C. D. (2002). Predictions of understorey light conditions in northern hardwood forests following parameterization, sensitivity analysis, and tests of the SORTIE light model. Forest Ecology and Management, 165, 235–248. https://doi.org/10.1016/S0378-1127(01)00621-1 Beaudet, M., Messier, C., & Leduc, A. (2004). Understorey light profiles in temperate deciduous forests: Recovery process following selection cutting. Journal of Ecology, 92(2), 328– 338. https://doi.org/10.1111/j.0022-0477.2004.00869.x Bennett, J. A., Maherali, H., Reinhart, K. O., Lekberg, Y., Hart, M. M., & Klironomos, J. (2017). Plant-soil feedbacks and mycorrhizal type influence temperate forest population dynamics. In Science (Vol. 355). https://www.science.org Berkowitz, A. R., Canham, C. D., & Kelly, V. R. (1995). Competition vs. facilitation of tree seedling growth and survival in early successional communities. Ecology, 76(4), 1156– 1168. https://doi.org/10.2307/1940923 Bevilacqua, E., Nyland, R. D., Namestnik, T. S., & Allen, D. C. (2021). Growth of sugar maple (Acer saccharum marsh.) after ice storm damage and forest tent caterpillar (Malacosoma disstria hubner) defoliation. Forests, 12(5). https://doi.org/10.3390/f12050620 Beyer, D., Rudolph, B. A., Kintigh, K., Albright, C., Swanson, K., Smith, L. M., Begalle, D., & Doepker, R. (2010). Habitat and behavior of wintering deer in northern Michigan: a glossary of terms and associated background information. Michigan Department of Natural Resources and Environment. Wildlife Division Report No. 3520, November., 22. Bjorkbom, J. C. (1979). Seed production and advance regeneration in Allegheny hardwood forests. Res. Pap. NE-435. Broomall, PA: U.S. Department of Agriculture, Forest Service, Northeastern Forest Experiment Station. 10p. Bohlen, P. J., Scheu, S., Hale, C. M., McLean, M. A., Migge, S., Groffman, P. M., & Parkinson, D. (2004). Invasive earthworms in temperate forests. Frontiers in Ecology and Evolution, 2(8), 427–435. 152 Bohn, K. K., & Nyland, R. D. (2003). Forecasting development of understory American beech after partial cutting in uneven-aged northern hardwood stands. Forest Ecology and Management, 180, 453–461. https://doi.org/10.1016/S0378-1127(02)00614-X Bolton, N. W., & D’Amato, A. W. (2011). Regeneration responses to gap size and coarse woody debris within natural disturbance-based silvicultural systems in northeastern Minnesota, USA. Forest Ecology and Management, 262(7), 1215–1222. https://doi.org/10.1016/j.foreco.2011.06.019 Bradshaw, L., & Waller, D. M. (2016). Impacts of white-tailed deer on regional patterns of forest tree recruitment. Forest Ecology and Management, 375, 1–11. https://doi.org/10.1016/j.foreco.2016.05.019 Braga, S. R., de Oliveira, M. L. R., & Gorgens, E. B. (2021). forestmangr: Forest Mensuration and Management. R package version 0.9.3. https://CRAN.R- project.org/package=forestmangr. Braun, E. L. (1950). Deciduous forests of eastern North America. Philadelphia, PA: Blakiston. 596p. Brinkman, T. J., Person, D. K., Smith, W., Chapin, F. S., McCoy, K., Leonawicz, M., & Hundertmark, K. J. (2013). Using DNA to test the utility of pellet-group counts as an index of deer counts. Wildlife Society Bulletin, 37(2), 444–450. https://doi.org/10.1002/wsb.270 Brockerhoff, E. G., Barbaro, L., Castagneyrol, B., Forrester, D. I., Gardiner, B., González- Olabarria, J. R., Lyver, P. O. B., Meurisse, N., Oxbrough, A., Taki, H., Thompson, I. D., van der Plas, F., & Jactel, H. (2017). Forest biodiversity, ecosystem functioning and the provision of ecosystem services. Biodiversity and Conservation, 26(13), 3005–3035. https://doi.org/10.1007/s10531-017-1453-2 Burger, T. L., & Kotar, J. (2003). A guide to forest communities and habitat types of Michigan. University of Wisconsin-Madison. Burns, R. M., & Honkala, B. H. (1990). Silvics of North America: 2. Hardwoods. Agriculture Handbook 654. United States Department of Agriculture, Forest Service. Byun, K., & Hamlet, A. F. (2018). Projected changes in future climate over the Midwest and Great Lakes region using downscaled CMIP5 ensembles. International Journal of Climatology, 38, e531–e553. https://doi.org/10.1002/joc.5388 Campbell, T. A., Laseter, B. R., Ford, W. M., & Miller, K. V. (2004). Movements of female white-tailed deer (Odocoileus virginianus) in relation to timber harvests in the central Appalachians. Forest Ecology and Management, 199(2–3), 371–378. https://doi.org/10.1016/j.foreco.2004.05.051 Canham, C. D. (1985). Suppression and release during canopy recruitment in Acer saccharum. Torrey Botanical Club, 112(2), 134–145. 153 Canham, C. D., & Burbank, D. H. (1994). Causes and consequences of resource heterogeneity in forests: interspecific variation in light transmission by canopy trees. Canadian Journal of Forest Research, 24, 337–349. Caspersen, J. P., & Kobe, R. K. (2001). Interspecific variation in sapling mortality in relation to growth and soil moisture. Oikos, 92(1), 160–168. Caspersen, J. P., & Saprunoff, M. (2005). Seedling recruitment in a northern temperate forest: The relative importance of supply and establishment limitation. Canadian Journal of Forest Research, 35(4), 978–989. https://doi.org/10.1139/x05-024 Chao, A., Gotelli, N. J., Hsieh, T. C., Sander, E. L., Ma, K. H., Colwell, R. K., & Ellison, A. M. (2014). Rarefaction and extrapolation with Hill numbers: A framework for sampling and estimation in species diversity studies. Ecological Monographs, 84(1), 45–67. https://doi.org/10.1890/13-0133.1 Chapman, D. G. (1961). Statistical problems in dynamics of exploited fisheries populations. Fourth Berkeley Symposium, 153–168. Chapman, D., Purse, B. v., Roy, H. E., & Bullock, J. M. (2017). Global trade networks determine the distribution of invasive non-native species. Global Ecology and Biogeography, 26(8), 907–917. https://doi.org/10.1111/geb.12599 Cheng, S., Widden, P., & Messier, C. (2005). Light and tree size influence belowground development in yellow birch and sugar maple. Plant and Soil, 270(1), 321–330. https://doi.org/10.1007/s11104-004-1726-x Cleavitt, N. L., Battles, J. J., Johnson, C. E., & Fahey, T. J. (2018). Long-term decline of sugar maple following forest harvest, Hubbard Brook Experimental Forest, New Hampshire. Canadian Journal of Forest Research, 48(1), 23–31. https://doi.org/10.1139/cjfr-2017- 0233 Cleavitt, N. L., & Fahey, T. J. (2017). Seed production of sugar maple and American beech in northern hardwood forests, New Hampshire, USA. Canadian Journal of Forest Research, 47, 985–990. https://doi.org/10.1139/cjfr-2017-0096 Cole, W. G., & Lorimer, C. G. (2005). Probabilities of small-gap capture by sugar maple saplings based on height and crown growth data from felled trees. Canadian Journal of Forest Research, 35, 643–655. https://doi.org/10.1139/X04-210 Collin, A., Messier, C., Kembel, S. W., & Bélanger, N. (2017). Low light availability associated with American beech is the main factor for reduced sugar maple seedling survival and growth rates in a hardwood forest of Southern Quebec. Forests, 8(11), 1–13. https://doi.org/10.3390/f8110413 Côté, S. D., Rooney, T. P., Tremblay, J., Dussault, C., & Waller, D. M. (2004). Ecological impacts of deer overabundance. Annual Review of Ecology, Evolution, and Systematics, 35, 113–147. https://doi.org/10.2307/annurev.ecolsys.35.021103.30000006 154 Crow, T. R., Buckley, D. S., Nauertz, E. A., & Zasada, J. C. (2002). Effects of management on the composition and structure of northern hardwood forests in upper Michigan. Forest Science, 48(1), 129–145. https://doi.org/10.1093/forestscience/48.1.129 Curtis, J. T., & McIntosh, R. P. (1951). An upland forest continuum in the prairie-forest border region of Wisconsin. Ecology, 32(3), 476–496. https://doi.org/10.2307/1931725 Curtis, R., & Rushmore, F. (1958). Some effects of stand density and deer browsing on reproduction in an Adirondack hardwood stand. Journal of Forestry, 56(2), 116–121. https://doi.org/10.1093/jof/56.2.116 Dale, V. H., Joyce, L. A., McNulty, S., & Neilson, R. P. (2000). The interplay between climate change, forests, and disturbances. Science of the Total Environment, 262(3), 201–204. https://doi.org/10.1016/S0048-9697(00)00522-2 Dale, V. H., Joyce, L. A., McNulty, St., Neilson, R. P., Ayres, M. P., Flannigan, M. D., Hanson, P. J., Irland, L. C., Lugo, A. E., Peterson, C. J., Simberloff, D., Swanson, F. J., Stocks, B. J., & Wotton, B. M. (2001). Climate change and forest disturbances. BioScience, 51(9), 723–734. Danyagri, G., Baral, S. K., & Pelletier, G. (2019). Effects of disturbance and site factors on sapling dynamics and species diversity in northern hardwood stands. Forest Ecology and Management, 444(May), 225–234. https://doi.org/10.1016/j.foreco.2019.04.041 Dawe, K. L., Bayne, E. M., & Boutin, S. (2014). Influence of climate and human land use on the distribution of white-tailed deer (Odocoileus virginianus) in the western boreal forest. Canadian Journal of Zoology, 92, 353–363. Dey, D. C. (2014). Sustaining oak forests in eastern North America: Regeneration and recruitment, the pillars of sustainability. Forest Science, 60(5), 926–942. https://doi.org/10.5849/forsci.13-114 Dey, D. C., Dwyer, J., & Wiedenbeck, J. (2017). Relationship between tree value, diameter, and age in high-quality sugar maple (Acer saccharum) on the Menominee Reservation, Wisconsin. Journal of Forestry, 115(5), 397–405. Diamond, J. (1986). Overview: Laboratory experiments, field experiments, and natural experiments. In J. Diamond & T. J. Case (Eds.), Community Ecology (pp. 3–22). Harper & Row. Dickmann, D. I., & Leefers, L. A. (2016). The Forests of Michigan (2nd ed.). University of Michigan Press. Donoso, P. J., Nyland, R. D., & Zhang, L. (2000). Growth of saplings after selection cutting in northern hardwoods. Northern Journal of Applied Forestry, 17(4), 149–152. 155 Donovan, G. (2005). Chronic regeneration failure in northern hardwood stands: A liability to certified forest landowners. Proceedings of the Forests & Wildlife-Striving for Balance Michigan Society of American Foresters Conference, 125–130. Dukes, J. S., Pontius, J., Orwig, D., Garnas, J. R., Rodgers, V. L., Brazee, N., Cooke, B., Theoharides, K. A., Stange, E. E., Harrington, R., Ehrenfeld, J., Gurevitch, J., Lerdau, M., Stinson, K., Wick, R., & Ayres, M. (2009). Responses of insect pests, pathogens, and invasive plant species to climate change in the forests of northeastern North America: What can we predict? Canadian Journal of Forest Research, 39, 231–248. Duncan, R. P. (1989). An evaluation of errors in tree age estimates based on increment cores in kahikatea (Dacrycarpus dacrydioides). New Zealand Natural Sciences, 16(January), 31– 37. Duval, R. P., McConnell, T. E., & Hix, D. M. (2014). Annual change in Ohio hardwood stumpage prices, 1960 to 2011. Forest Products Journal, 64(1/2), 19–25. Elenitsky, L. M., Walters, M. B., & Farinosi, E. J. (2020). Tree regeneration structure following beech bark disease-motivated harvests: Factors associated with patterns and management implications. Forests, 11(2), 180. https://doi.org/10.3390/f11020180 Elzhov, T. v., Mullen, K. M., Spiess, A.-N., & Bolker, B. (2016). minpack.lm: R Interface to the Levenberg-Marquardt Nonlinear Least-Squares Algorithm Found in MINPACK, Plus Support for Bounds. R package version 1.2-1. Eyre, F. H., & Zillgitt, W. M. (1953). Partial cuttings in northern hardwoods of the Lake States: twenty-year experimental results. Technical Bulletin LS-1076, 124. https://www.nrs.fs.fed.us/pubs/4626 Fassnacht, K. S., & Gower, S. T. (1997). Interrelationships among the edaphic and stand characteristics, leaf area index, and aboveground net primary production of upland forest ecosystems in north central Wisconsin. Canadian Journal of Forest Research, 27, 1058– 1067. https://doi.org/10.1139/x97-058 Fei, S., & Steiner, K. C. (2007). Evidence for increasing red maple abundance in the eastern United States. Forest Science, 53(4), 473–477. https://doi.org/10.1093/forestscience/53.4.473 Finley, A. O., & Banerjee, S. (2020). Bayesian spatially varying coefficient models in the spBayes R package. Environmental Modelling and Software, 125, 104608. https://doi.org/10.1016/j.envsoft.2019.104608 Finley, A. O., Banerjee, S., & Gelfand, A. E. (2015). spBayes for large univariate and multivariate point-referenced spatio-temporal data models. Journal of Statistical Software, 63(13), 1–28. https://doi.org/10.18637/jss.v063.i13 156 Flannigan, M. D., Stocks, B. J., & Wotton, B. M. (2000). Climate change and forest fires. Science of The Total Environment, 262, 221–229. https://doi.org/10.1016/S0048- 9697(00)00524-6 Forrester, J. A., Lorimer, C. G., Dyer, J. H., Gower, S. T., & Mladenoff, D. J. (2014). Response of tree regeneration to experimental gap creation and deer herbivory in north temperate forests. Forest Ecology and Management, 329, 137–147. https://doi.org/10.1016/j.foreco.2014.06.025 Forrester, J. A., Mcgee, G. G., & Mitchell, M. J. (2003). Effects of beech bark disease on aboveground biomass and species composition in a mature northern hardwood forest. The Journal of the Torrey Botanical Society, 130(2), 70–78. Forsyth, D. M., Barker, R. J., Moriss, G., & Scroggie, M. P. (2007). Modeling the relationship between fecal pellet indices and deer density. The Journal of Wildlife Management, 71(3), 964–970. https://doi.org/10.2193/2005-695 Fox, J., & Weisberg, S. (2011). An {R} Companion to Applied Regression, Second Edition. Thousand Oaks CA: Sage. URL: http://socserv.socsci.mcmaster.ca/jfox/Books/Companion. Frelich, L. E., & Lorimer, C. G. (1991). Natural disturbance regimes in hemlock-hardwood forests of the upper Great Lakes region. Ecological Monographs, 61(2), 145–164. https://doi.org/10.2307/1943005 Garrett, P., & Graber, R. (1995). Sugar maple seed production in northern New Hampshire. USDA Forest Service Northeastern Forest Experiment Station, 697, 1–6. Gelman, A., & Hill, J. (2007). Data analysis using regression and multilevel/hierarchical models. Cambridge University Press. Global Forest Resources Assessment 2015. (2015). In Food and Agriculture Organization of the United Nations. https://doi.org/10.1002/2014GB005021 Goldblum, D., & Rigg, L. S. (2002). Age structure and regeneration dynamics of sugar maple at the deciduous/boreal forest ecotone, Ontario, Canada. Physical Geography, 23(2), 115– 129. https://doi.org/10.2747/0272-3646.23.2.115 Gravel, D., Beaudet, M., & Messier, C. (2011). Sapling age structure and growth series reveal a shift in recruitment dynamics of sugar maple and American beech over the last 40 years. Canadian Journal of Forest Research, 41(4), 873–880. https://doi.org/10.1139/x10-242 Grayson, S. F., Buckley, D. S., Henning, J. G., Schweitzer, C. J., Gottschalk, K. W., & Loftis, D. L. (2012). Understory light regimes following silvicultural treatments in central hardwood forests in Kentucky, USA. Forest Ecology and Management, 279, 66–76. https://doi.org/10.1016/j.foreco.2012.05.017 157 Grisez, T. J. (1960). Slash helps protect seedlings from deer browsing. Journal of Forestry, 58(5), 385–387. Hagge, J., Müller, J., Bässler, C., Biebl, S. S., Brandl, R., Drexler, M., Gruppe, A., Hotes, S., Hothorn, T., Langhammer, P., Stark, H., Wirtz, R., Zimmerer, V., & Mysterud, A. (2019). Deadwood retention in forests lowers short-term browsing pressure on silver fir saplings by overabundant deer. Forest Ecology and Management, 451, 117531. https://doi.org/10.1016/j.foreco.2019.117531 Hale, C. M., Frelich, L. E., & Reich, P. B. (2006). Changes in hardwood forest understory plant communities in response to European earthworm invasions. Ecology, 87(7), 1637–1649. http://climate.umn.edui Hanberry, B. B., & Abrams, M. D. (2019). Does white-tailed deer density affect tree stocking in forests of the Eastern United States? Ecological Processes, 8(30). https://doi.org/10.1186/s13717-019-0185-5 Hane, E. N. (2003). Indirect effects of beech bark disease on sugar maple seedling survival. Canadian Journal of Forest Research, 33, 807–813. https://doi.org/10.1139/x03-008 Harmala, H. (2021). Understory dynamics across 62-years of a northern hardwood management gradient study (Issue September). Michigan Technological University. Henne, P. D., Hu, F. S., & Cleland, D. T. (2007). Lake-effect snow as the dominant control of mesic-forest distribution in Michigan, USA. Journal of Ecology, 95, 517–529. https://doi.org/10.1111/j.1365-2745.2007.01220.x Henry, C. R., Walters, M. B., Finley, A. O., Roloff, G. J., & Farinosi, E. J. (2021). Complex drivers of sugar maple (Acer saccharum) regeneration reveal challenges to long-term sustainability of managed northern hardwood forests. Forest Ecology and Management, 479(June 2020), 118541. https://doi.org/10.1016/j.foreco.2020.118541 Hett, J. M. (1971). A dynamic analysis of age in sugar maple seedlings. Ecology, 52(6), 1071– 1074. https://doi.org/10.2307/1933815 Holdsworth, A. R., Frelich, L. E., & Reich, P. B. (2007). Effects of earthworm invasion on plant species richness in northern hardwood forests. Conservation Biology, 21(4), 997–1008. https://doi.org/10.1111/j.1523-1739.2007.00740.x Holdsworth, A. R., Frelich, L. E., & Reich, P. B. (2008). Litter decomposition in earthworm- invaded northern hardwood forests: Role of invasion degree and litter chemistry. Ecoscience, 15(4), 536–544. https://doi.org/10.2980/15-4-3151 Holter, J. B., Urban Jr., W. E., Hayes, H. H., Silver, H., & Skutt, H. R. (1975). Ambient temperature effects on physiological traits of white-tailed deer. Canadian Journal of Zoology, 53(6), 679–685. 158 Horsley, S. B., Stout, S. L., & DeCalesta, D. S. (2003). White-tailed deer impact on the vegetation dynamics of a northern hardwood forest. Ecological Applications, 13(1), 98– 118. Hsieh, T. C., Ma, K. H., & Chao, A. (2016). iNEXT: an R package for rarefaction and extrapolation of species diversity (Hill numbers). Methods in Ecology and Evolution, 7(12), 1451–1456. https://doi.org/10.1111/2041-210X.12613 Hughes, J. W., & Fahey, T. J. (1991). Colonization dynamics of herbs and shrubs in a disturbed northern hardwood forest. Source: Journal of Ecology, 79(3), 605–616. Hupperts, S. F., Webster, C. R., Froese, R. E., & Dickinson, Y. L. (2020). Seedling and sapling recruitment following novel silvicultural treatments in Great Lakes northern hardwoods. Forest Ecology and Management, 462(December 2019), 117983. https://doi.org/10.1016/j.foreco.2020.117983 Iverson, L. R., Schwartz, M. W., & Prasad, A. M. (2004). How fast and far might tree species migrate in the eastern United States due to climate change? Global Ecology and Biogeography, 13(3), 209–219. https://doi.org/10.1111/j.1466-822X.2004.00093.x James, G., Witten, D., Hastie, T., & Tibshirani, R. (Eds.). (2013). An introduction to statistical learning: with applications in R. Springer. Jenkins, J. (1997). Hardwood regeneration failure in the Adirondacks: Preliminary studies of incidence and severity. Wildlife Conservation Society 9(9), 70p. Kain, M., Battaglia, L., Royo, A., & Carson, W. P. (2011). Over-browsing in Pennsylvania creates a depauperate forest dominated by an understory tree: Results from a 60-year-old deer exclosure. Journal of the Torrey Botanical Society, 138(3), 322–326. https://doi.org/10.3159/TORREY-D-11-00018.1 Kassambara, A. (2020). rstatix: Pipe-friendly framework for basic statistical tests. R package version 0.5.0. https://CRAN.R-project.org/package=rstatix. Kern, C. C., D’Amato, A. W., & Strong, T. F. (2013). Diversifying the composition and structure of managed, late-successional forests with harvest gaps: What is the optimal gap size? Forest Ecology and Management, 304, 110–120. https://doi.org/10.1016/j.foreco.2013.04.029 Kern, C., Erdmann, G., Kenefic, L., Palik, B., & Strong, T. (2014). Development of the selection system in northern hardwood forests of the Lake States: An 80-year silviculture research legacy. In USDA Forest Service Experimental Forests and Ranges: Research for the Long Term (pp. 201–223). https://doi.org/10.1007/978-1-4614-1818-4 Kern, C. C., Burton, J. I., Raymond, P., D’Amato, A. W., Keeton, W. S., Royo, A. A., Walters, M. B., Webster, C. R., & Willis, J. L. (2017). Challenges facing gap-based silviculture and possible solutions for mesic northern forests in North America. Forestry, 90, 4–17. https://doi.org/10.1093/forestry/cpw024 159 Klooster, W. S., Herms, D. A., Knight, K. S., Herms, C. P., McCullough, D. G., Smith, A., Gandhi, K. J. K., & Cardina, J. (2014). Ash (Fraxinus spp.) mortality, regeneration, and seed bank dynamics in mixed hardwood forests following invasion by emerald ash borer (Agrilus planipennis). Biological Invasions, 16, 859–873. https://doi.org/10.1007/s10530- 013-0543-7 Knapp, S. P., Kern, C. C., & Webster, C. R. (2021). Harvested opening size affects cohort development and failures in a second-growth northern hardwood forest. Forest Ecology and Management, 482(September 2020), 118804. https://doi.org/10.1016/j.foreco.2020.118804 Kobe, R. K. (1996). Intraspecific variation in sapling mortality and growth predicts geographic variation in forest composition. Ecological Monographs, 66(2), 181–201. Kobe, R. K. (2006). Sapling growth as a function of light and landscape-level variation in soil water and foliar nitrogen in northern Michigan. Oecologia, 147, 119–133. Koorem, K., Tulva, I., Davison, J., Jairus, T., Öpik, M., Vasar, M., Zobel, M., & Moora, M. (2017). Arbuscular mycorrhizal fungal communities in forest plant roots are simultaneously shaped by host characteristics and canopy-mediated light availability. Plant and Soil, 410(1–2), 259–271. https://doi.org/10.1007/s11104-016-3004-0 Leak, W. B., & Wilson Jr., R. W. (1958). Regeneration after cutting of old-growth northern hardwoods in New Hampshire. USDA Forest Service, Northeast Forest Experimental Station, Station Paper, 6, 8. Leak, W. B. (2006). Sixty years of change in the sapling understories of northern hardwood selection stands in New Hampshire. Northern Journal of Applied Forestry, 23(4), 301– 303. Lenth, R. (2020). emmeans: Estimated Marginal Means, aka Least-Squares Means. R package version 1.5.2-1. https://CRAN.R-project.org/package=emmeans. Levenberg, K. (1944). A method for the solution of certain non-linear problems in least squares. The Quarterly of Applied Mathematics, 2, 164–169. Liebhold, A. M., Mccullough, D. G., Blackburn, L. M., Frankel, S. J., von Holle, B., & Aukema, J. E. (2013). A highly aggregated geographical distribution of forest pest invasions in the USA. Diversity and Distributions, 19(9), 1208–1216. https://doi.org/10.1111/ddi.12112 Linehan, P. E., & Jacobson, M. G. (2005). Forecasting hardwood stumpage price trends in Pennsylvania. Forest Products Journal, 55(12), 47–52. Lorimer, C. G., Dahir, S. E., & Singer, M. T. (1999). Frequency of partial and missing rings in Acer saccharum in relation to canopy position and growth rate. Plant Ecology, 143, 189– 202. https://doi.org/10.1023/A:1009847819158 160 Lüdecke; Makowski; Waggoner; Patil. (2020). Assessment of Regression Models Performance. CRAN. Available from https://easystats.github.io/performance/. Malis, F., Kopecky, M., Petrik, P., Vladovic, J., Merganic, J., & Vida, T. (2016). Life stage, not climate change, explains observed tree range shifts. Global Change Biology, 22, 1904– 1914. https://doi.org/10.1111/gcb.13210 Mangiafico, S. (2021). rcompanion: Functions to support extension education program evaluation. https://CRAN.R-project.org/package=rcompanion. Marks, P. L., & Gardescu, S. (1998). A case study of sugar maple (Acer saccharum) as a forest seedling bank species. Journal of the Torrey Botanical Society, 125(4), 287–296. https://doi.org/10.2307/2997242 Marquardt, D. W. (1963). An algorithm for least-squares estimation of nonlinear parameters. Journal of the Society for Industrial and Applied Mathematics, 11(2), 431–441. Marquaridt, D. W. (1970). Generalized inverses, ridge regression, biased linear estimation, and nonlinear estimation. Technometrics, 12(3), 591–612. https://doi.org/10.1080/00401706.1970.10488699 Marquis, D. A. (1967). Clearcutting in Northern Hardwoods: Results after 30 years. Research Paper NE-85. Upper Darby, PA: U.S.D.A., Forest Service, Northern Research Station, 13. Marx, L., & Walters, M. B. (2008). Survival of tree seedlings on different species of decaying wood maintains tree distribution in Michigan hemlock-hardwood forests. Journal of Ecology, 96(3), 505–513. https://doi.org/10.1111/j.1365-2745.2008.01360.x Matonis, M. S., Walters, M. B., & Millington, J. D. A. (2011a). Gap-, stand-, and landscape- scale factors contribute to poor sugar maple regeneration after timber harvest. Forest Ecology and Management, 262(2), 286–298. https://doi.org/10.1016/j.foreco.2011.03.034 McEuen, A. B., & Curran, L. M. (2004). Seed dispersal and recruitment limitation across spatial scales in temperate forest fragments. Ecology, 85(2), 507–518. McWilliams, W. H., Stout, S. L., Bowersox, T. W., & McCormick, L. H. (1995). Adequacy of advance tree-seedling regeneration in Pennsylvania’s forests. Northern Journal of Applied Forestry, 12(4), 187–191. https://doi.org/10.1093/njaf/12.4.187 MDNR. (2016). Michigan deer management plan. Michigan Department of Natural Resources, Wildlife Division Report No. 3626. Metzger, F. T., & Tubbs, C. H. (1971). The influence of cutting method on regeneration of second-growth northern hardwoods. Journal of Forestry, 69(9), 559–564. 161 Millar, C. I., Stephenson, N. L., & Stephens, S. L. (2007). Climate change and forest of the future: Managing in the face of uncertainty. Ecological Applications, 17(8), 2145–2151. https://doi.org/http://dx.doi.org/10.1890/06-1715.1 Miller, K. M., & McGill, B. J. (2019). Compounding human stressors cause major regeneration debt in over half of eastern US forests. Journal of Applied Ecology, 56(6), 1355–1366. https://doi.org/10.1111/1365-2664.13375 Miller, R. O. (2004). Regeneration in a heavily browsed northern hardwood stand twelve years after scarification and fencing. Millington, J. D. A., Walters, M. B., Matonis, M. S., & Liu, J. (2010). Effects of local and regional landscape characteristics on wildlife distribution across managed forests. Forest Ecology and Management, 259, 1102–1110. https://doi.org/10.1016/j.foreco.2009.12.020 Millington, J. D. A., Walters, M. B., Matonis, M. S., Laurent, E. J., Hall, K. R., & Liu, J. (2011). Combined long-term effects of variable tree regeneration and timber management on forest songbirds and timber production. Forest Ecology and Management, 262(5), 718– 729. https://doi.org/10.1016/j.foreco.2011.05.002 Minor, D. M., & Kobe, R. K. (2017). Masting synchrony in northern hardwood forests: super- producers govern population fruit production. Journal of Ecology, 105(4), 987–998. https://doi.org/10.1111/1365-2745.12729 Moen, A. N. (1976). Energy conservation by white-tailed deer in the winter. Ecology, 57(1), 192–198. Moran, P. A. P. (1950). Notes on Continuous Stochastic Phenomena. Biometrika, 37, 17–23. National Operational Hydrologic Remote Sensing Center. (2004). Snow Data Assimilation System (SNODAS) Data Products at NSIDC, Version 1. Daily Snow Depth. NSIDC: National Snow and Ice Data Center. https://doi.org/https://doi.org/10.7265/N5TB14TC Neuendorff, J. K., Nagel, L. M., Webster, C. R., & Janowiak, M. K. (2007). Stand structure and composition in a northern hardwood forest after 40 years of single-tree selection. Northern Journal of Applied Forestry, 24(3), 197–202. Neuenkamp, L., Zobel, M., Koorem, K., Jairus, T., Davison, J., Öpik, M., Vasar, M., & Moora, M. (2021). Light availability and light demand of plants shape the arbuscular mycorrhizal fungal communities in their roots. Ecology Letters, 24(3), 426–437). https://doi.org/10.1111/ele.13656 Neumann, D. (2015). Silvics and Management Guidance Manual (Vol. 4111). Michigan Department of Natural Resources. Niese, J. N., & Strong, T. F. (1992). Economic and tree diversity trade-offs in managed northern hardwoods. Canadian Journal of Forest Research, 22, 1807–1813. 162 Nuckolls, A. E., Wurzburger, N., Ford, C. R., Hendrick, R. L., Vose, J. M., & Kloeppel, B. D. (2009). Hemlock declines rapidly with hemlock woolly adelgid infestation: Impacts on the carbon cycle of southern Appalachian forests. Ecosystems, 12, 179–190. https://doi.org/10.1007/s10021-008-9215-3 Nuttle, T., Royo, A. A., Adams, M. B., & Carson, W. P. (2013). Historic disturbance regimes promote tree diversity only under low browsing regimes in eastern deciduous forest. Ecological Monographs, 83(1), 3–17. https://doi.org/10.1890/11-2263.1 Nyland, R. D. (1992). Exploitation and Greed in Eastern Hardwood Forests. Journal of Forestry, January(1), 33–37. https://doi.org/10.5849/jof.12-060 Nyland, R. D. (2002). Silviculture: Concept and applications (2nd Edition). McGraw-Hill. Nyland, R. D., Nystrom, L., Kiernan, D. H., & Bevilacqua, E. (2019). Assessing small-stem density in northern hardwood selection system stands. Canadian Journal of Forest Research, 49, 237–245. https://doi.org/10.1139/cjfr-2018-0192 OMNRF. (2015). Forest management guide to silviculture in the Great Lakes-St. Lawrence and boreal forests of Ontario. Queens Printer for Ontario. Oswalt, S. N., Smith, W. B., Miles, P. D., & Pugh, S. A. (2014). Forest resources of the United States, 2012: A technical document supporting the Forest Service update of the 2010 RPA Assessment. In USDA General Technical Report WO-91. Ozoga, J. J., & Gysel, L. W. (1972). Response of white-tailed deer to winter weather. The Journal of Wildlife Management, 36(3), 892–896. Panzavolta, T., Bracalini, M., Benigno, A., & Moricca, S. (2021). Alien Invasive Pathogens and Pests Harming Trees, Forests, and Plantations: Pathways, Global Consequences and Management. Forests, 12(10), 1364. https://doi.org/10.3390/f12101364 Parker, G. R., & Leopold, D. J. (1983). Replacement of Ulmus americana L. in a mature east- central Indiana woods. Bulletin of the Torrey Botanical Club, 110(4), 482–488. Pellerin, M., Saïd, S., Richard, E., Hamann, J. L., Dubois-Coli, C., & Hum, P. (2010). Impact of deer on temperate forest vegetation and woody debris as protection of forest regeneration against browsing. Forest Ecology and Management, 260, 429–437. https://doi.org/10.1016/j.foreco.2010.04.031 Plummer, M. (2003). JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. Proceedings of the 3rd International Workshop on Distributed Statistical Computing, 1–10. https://doi.org/10.1002/ana.1067 Plummer, M., Best, N., Cowles, K., & Vines, K. (2006). CODA: Convergence diagnosis and output analysis for MCMC. R News, 6(1), 7–11. 163 Plummer, M. (2018). rjags: Bayesian graphical models using MCMC. R package version 4-8. https://CRAN.R-project.org/package=rjags. Polson, N. G., Scott, J. G., & Windle, J. (2013). Bayesian inference for logistic models using Pólya-Gamma latent variables. Journal of the American Statistical Association, 108(504), 1339–1349. https://doi.org/10.1080/01621459.2013.829001 Pond, N. C., Froese, R. E., & Nagel, L. M. (2014). Sustainability of the selection system in northern hardwood forests. Forest Science, 60(2), 374–381. Pond, N. C., & Froese, R. E. (2015). Interpreting stand structure through diameter distributions. Forest Science, 61(3), 429–437. https://doi.org/10.5849/forsci.14-056 Powers, M. D., & Nagel, L. M. (2009). Pennsylvania sedge cover, forest management and deer density influence tree regeneration dynamics in a northern hardwood forest. Forestry, 82(3), 241–254. https://doi.org/10.1093/forestry/cpp003 Poznanovic, S. K., Webster, C. R., & Bump, J. K. (2013). Maintaining mid-tolerant tree species with uneven-aged forest management: 9-year results from a novel group-selection experiment. Forestry, 86(5), 555–567. https://doi.org/10.1093/forestry/cpt025 R Core Team. (2018). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. Ramirez, J. I., Jansen, P. A., & Poorter, L. (2018). Effects of wild ungulates on the regeneration, structure and functioning of temperate forests: A semi-quantitative review. Forest Ecology and Management, 424(December 2017), 406–419. https://doi.org/10.1016/j.foreco.2018.05.016 Ramsfield, T. D., Bentz, B. J., Faccoli, M., Jactel, H., & Brockerhoff, E. G. (2016). Forest health in a changing world: Effects of globalization and climate change on forest insect and pathogen impacts. Forestry, 89(3), 245–252. https://doi.org/10.1093/forestry/cpw018 Randall, J. A., & Walters, M. B. (2011). Deer density effects on vegetation in aspen forest understories over site productivity and stand age gradients. Forest Ecology and Management, 261, 408–415. https://doi.org/10.1016/j.foreco.2010.10.026 Randall, J. A., & Walters, M. B. (2019). Competition and tolerance of low soil water favor Carex dominance over establishing Acer seedlings in managed temperate mesic forests. Forest Ecology and Management, 449(July), 117481. https://doi.org/10.1016/j.foreco.2019.117481 Resner, K., Yoo, K., Sebestyen, S. D., Aufdenkampe, A., Hale, C., Lyttle, A., & Blum, A. (2015). Invasive earthworms deplete key soil inorganic nutrients (Ca, Mg, K, and P) in a northern hardwood forest. Ecosystems, 18(1), 89–102. https://doi.org/10.1007/s10021- 014-9814-0 164 Richards, F. J. (1959). A flexible growth function for empirical use. Journal of Experimental Botany, 10, 290–300. https://academic.oup.com/jxb/article/10/2/290/528209 Rogers, B. M., Jantz, P., Goetz, S. J., & Hole, W. (2017). Vulnerability of eastern US tree species to climate change. Global Change Biology, 23, 3302–3320. https://doi.org/10.1111/gcb.13585 Rogers, N. S., D’amato, A. W., & Leak, W. B. (2021). Long-term evolution of composition and structure after repeated group selection over eight decades. Canadian Journal of Forest Research, 51(7), 1080–1091. https://doi.org/10.1139/cjfr-2020-0339 Royo, A. A., & Carson, W. P. (2006). On the formation of dense understory layers in forests worldwide: Consequences and implications for forest dynamics, biodiversity, and succession. Canadian Journal of Forest Research, 36(6), 1345–1362. https://doi.org/10.1139/X06-025 Sage, R. W., Porter, W. F., & Underwood, H. B. (2003). Windows of opportunity: White-tailed deer and the dynamics of northern hardwood forests of the northeastern US. Journal for Nature Conservation, 10, 213–220. https://doi.org/10.1078/1617-1381-00021 Sakamoto, Y., Ishiguro, M., & Kitagawa, G. (1986). Akaike Information Criterion Statistics. D. Reidel Publishing Company. Saunders, M. R., & Puettmann, K. J. (1999). Use of vegetational characteristics and browsing patterns to predict deer damage in eastern white pine (Pinus strobus) plantations. Northern Journal of Applied Forestry, 16(2), 96–102. https://doi.org/10.1093/njaf/16.2.96 Schomaker, M. E., Zarnoch, S. J., Bechtold, W. A., Latelle, D. J., Burkman, W. G., & Cox, S. M. (2007). Crown-Condition Classification: A Guide to Data Collection and Analysis. Gen Tech. Rep. SRS-102. Asheville, NC: U.S. Department of Agriculture, Forest Service, Southern Research Station. 78p. Schulte, L. A., Mladenoff, D. J., Crow, T. R., Merrick, L. C., & Cleland, D. T. (2007). Homogenization of northern U.S. Great Lakes forests due to land use. Landscape Ecology, 22, 1089–1103. https://doi.org/10.1007/s10980-007-9095-5 Schwarz, P. A., Fahey, T. J., & Mcculloch, C. E. (2003). Factors controlling spatial variation of tree species abundance in a forested landscape. Ecology, 84(7), 1862–1878. Schwinning, S., & Weiner, J. (1998). Mechanisms determining the degree of size asymmetry in competition among plants. Oecologia, 113(4), 447–455. https://doi.org/10.1007/s004420050397 Sendall, K. M., Lusk, C. H., & Reich, P. B. (2015). Becoming less tolerant with age: Sugar maple, shade, and ontogeny. Oecologia, 179(4), 1011–1021. https://doi.org/10.1007/s00442-015-3428-x 165 Shi, H., Laurent, E. J., LeBouton, J., Racevskis, L., Hall, K. R., Donovan, M., Doepker, R. V., Walters, M. B., Lupi, F., & Liu, J. (2006). Local spatial modeling of white-tailed deer distribution. Ecological Modelling, 190, 171–189. https://doi.org/10.1016/j.ecolmodel.2005.04.007 Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & van der Linde, A. (2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society. Series B: Statistical Methodology, 64(4), 583–616. https://doi.org/10.1111/1467-9868.00353 Telfer, E. S. (1970). Winter habitat selection by moose and white-tailed deer. The Journal of Wildlife Management, 34(3), 553–559. Thompson, I., Mackey, B., McNulty, S., & Mosseler, A. (2009). Forest resilience, biodiversity, and climate change. A synthesis of the biodiversity/resilience/stability relationship in forest ecosystems. Secretariat of the Convention on Biological Diversity, Montreal. Technical series. No 43. 67p. Tremblay, J. P., Huot, J., & Potvin, F. (2006). Divergent nonlinear responses of the boreal forest field layer along an experimental gradient of deer densities. Oecologia, 150(1), 78–88. https://doi.org/10.1007/s00442-006-0504-2 Tubbs, C. H. (1977a). Age and structure of a northern hardwood selection forest, 1929-1976. Journal of Forestry, 75(1), 22–24. Tubbs, C. H. (1977b). Manager's handbook for northern hardwoods in the north central states. U.S. Department of Agriculture Forest Service. Urbanek, R. E., Nielsen, C. K., Preuss, T. S., & Glowacki, G. A. (2012). Comparison of aerial surveys and pellet-based distance sampling methods for estimating deer density. Wildlife Society Bulletin, 36(1), 100–106. https://doi.org/10.1002/wsb.116 Van Deelen, T. R., Campa III, H., Hamady, M., & Haufler, J. B. (1998). Migration and seasonal range dynamics of deer using adjacent deeryards in northern Michigan. The Journal of Wildlife Management, 62(1), 205–213. Venables, W. N., & Ripley, B. D. (2002). Modern applied statistics with S. (Fourth Edition). Springer. Vickers, L. A., McWilliams, W. H., Knapp, B. O., D’Amato, A. W., Dey, D. C., Dickinson, Y. L., Kabrick, J. M., Kenefic, L. S., Kern, C. C., Larsen, D. R., Royo, A. A., Saunders, M. R., Shifley, S. R., & Westfall, J. A. (2019). Are current seedling demographics poised to regenerate northern US forests? Journal of Forestry, 117(6), 592–612. https://doi.org/10.1093/jofore/fvz046 Walters, M. B., & Reich, P. B. (1997). Growth of Acer saccharum seedlings in deeply shaded understories of northern Wisconsin: Effects of nitrogen and water availability. Canadian Journal of Forest Research, 27, 237–247. https://doi.org/10.1139/x96-178. 166 Walters, M. B., Farinosi, E. J., Willis, J. L., & Gottschalk, K. W. (2016). Managing for diversity: Harvest gap size drives complex light, vegetation, and deer herbivory impacts on tree seedlings. Ecosphere, 7(8), e01397. https://doi.org/10.1002/ecs2.1397 Walters, M. B., Farinosi, E. J., & Willis, J. L. (2020a). Deer browsing and shrub competition set sapling recruitment height and interact with light to shape recruitment niches for temperate forest tree species. Forest Ecology and Management, 467, 118134. https://doi.org/10.1016/j.foreco.2020.118134 Walters, M. B., Roloff, G. J., Henry, C. R., Hartman, J. P., Donovan, M. L., Farinosi, E. J., & Starking, M. D. (2020b). Rethinking northern hardwood forest management paradigms with silvicultural systems research: Research–management partnerships ensure relevance and application. Journal of Forestry, 118(3), 260–274. https://doi.org/10.1093/jofore/fvz071 Walters, M. B., Henry, C. R., Farinosi, E. J., Roloff, G. J., Donovan, M. L., & Hartman, J. P. (2022). Sapling stocking targets for multiple management goals in northern hardwood forests: How do stands measure up? Journal of Forestry, In press., 1–25. https://doi.org/10.1093/jofore/fvac002 Webster, C. R., Dickinson, Y. L., Burton, J. I., Frelich, L. E., Jenkins, M. A., Kern, C. C., Raymond, P., Saunders, M. R., Walters, M. B., & Willis, J. L. (2018). Promoting and maintaining diversity in contemporary hardwood forests: Confronting contemporary drivers of change and the loss of ecological memory. Forest Ecology and Management, 421, 98–108. https://doi.org/10.1016/j.foreco.2018.01.010 Wetzel, J. F., Wambaugh, J. R., & Peek, J. M. (1975). Appraisal of white-tailed deer winter habitats in northeastern Minnesota. The Journal of Wildlife Management, 39(1), 59–66. White, M. A., & Mladenoff, D. J. (1994). Old-growth forest landscape transitions from pre- European settlement to present. Landscape Ecology, 9(3), 191–205. https://doi.org/10.1007/BF00134747 White, M. A. (2012). Long-term effects of deer browsing: Composition, structure and productivity in a northeastern Minnesota old-growth forest. Forest Ecology and Management, 269, 222–228. https://doi.org/10.1016/j.foreco.2011.12.043 Whitney, G. G. (1987). An ecological history of the Great Lakes forest of Michigan. Journal of Ecology, 75(3), 667–684. Willis, J. L., Walters, M. B., & Gottschalk, K. W. (2015). Scarification and gap size have interacting effects on northern temperate seedling establishment. Forest Ecology and Management, 347, 237–246. https://doi.org/10.1016/j.foreco.2015.02.026 Willis, J. L., Walters, M., & Farinosi, E. (2016). Local seed source availability limits young seedling populations for some species more than other factors in northern hardwood forests. Forest Science, 62(4), 440–448. https://doi.org/10.5849/forsci.15-143 167 Yañez-Arenas, C., Martínez-Meyer, E., Mandujano, S., & Rojas-Soto, O. (2012). Modelling geographic patterns of population density of the white-tailed deer in central Mexico by implementing ecological niche theory. Oikos, 121, 2081–2089. https://doi.org/10.1111/j.1600-0706.2012.20350.x Zak, D. R., Pregitzer, K. S., & Host, G. E. (1986). Landscape variation in nitrogen mineralization and nitrification. Canadian Journal of Forest Research, 16(6), 1258–1263. https://doi.org/10.1139/x86-223 Zak, D. R., Host, G. E., & Pregitzer, K. S. (1989). Regional variability in nitrogen mineralization, nitrification, and overstory biomass in northern Lower Michigan. Canadian Journal of Forest Research, 19(12), 1521–1526. https://doi.org/10.1139/x89- 231 Zhu, K., Woodall, C. W., & Clark, J. S. (2012). Failure to migrate: Lack of tree range expansion in response to climate change. Global Change Biology, 18, 1042–1052. https://doi.org/10.1111/j.1365-2486.2011.02571.x 168