MODELING INDIVIDUAL VARIABILITY IN GROWTH AND ITS IMPORTANCE: AN APPLICATION FOR LAKE TROUT (SALVELINUS NAMAYCUSH) IN LAKE SUPERIOR By Elizabeth Stebbins A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Fisheries and Wildlife – Master of Science 2022 ABSTRACT MODELING INDIVIDUAL VARIABILITY IN GROWTH AND ITS IMPORTANCE: AN APPLICATION FOR LAKE TROUT (SALVELINUS NAMAYCUSH) IN LAKE SUPERIOR By Elizabeth Stebbins Correctly characterizing growth of fish within a population is a crucial component of fish biology and fishery management because, among other things, it informs population dynamics that affect management decisions. Size-at-age is a common metric of fish growth and is often measured at the population level with the assumption that, on average, all fish of a given age are a given size. Over time, several studies have shown that ignoring individual variability in growth can influence population parameter estimates and these inaccuracies can be propagated in population models that are used to calculate reference points for management. In the first chapter we develop a hierarchical, mixed-effects statistical growth model that measures individual variability in growth model parameters and partitions it into two sources. We fit this model to length-at-age data of lake trout (Salvelinus namaycush) from six populations in Lake Superior and show that individual-level variability exceeds population-level variability for this system, and persistent error contributes more to variability in length-at-age. In our second chapter, we simulate a population of fish and predict biological reference points, yield-per-recruit, and spawning stock biomass-per-recruit curves from the population using a ‘standard’ method that ignores individual variability and a ‘true’ method that accounts for size-selective mortality and its interaction with individual fish. We show that ignoring individual variability in these models results in overestimation of yield-per-recruit and the biological reference points F0.1 and FMAX. Further, spawning stock biomass-per-recruit is underestimated at low levels of fishing intensity and overestimated at high levels of fishing intensity when individual variability is ignored. This thesis is dedicated to my parents, my siblings, my family, and my dear friends. (And to the COVID-19 vaccine.) iii ACKNOWLEDGEMENTS This work is supported by the Great Lakes Fishery Commission Project 2019_BEN_44079. I would like to thank Steve Lenart, Jason Smith, and Shawn Sitar for providing context on modeling and management of the Great Lakes. iv TABLE OF CONTENTS LIST OF TABLES ........................................................................................................................ vii LIST OF FIGURES ..................................................................................................................... viii INTRODUCTION ...........................................................................................................................1 CHAPTER ONE: MODELING INDIVIDUAL VARIABILITY IN LAKE TROUT GROWTH ............ 5 ABSTRACT ................................................................................................................................5 INTRODUCTION ......................................................................................................................6 METHODS ...............................................................................................................................11 Data Collection ....................................................................................................................11 The von Bertalanffy Growth Model for an Individual Fish .................................................11 Hierarchical model structure ...............................................................................................13 Distribution of observed data ..............................................................................................14 Estimation ............................................................................................................................14 Quantification of Persistent vs. Transient Variation ...........................................................15 Alternative Models: Biphasic and Sex-Specific ...................................................................16 RESULTS .................................................................................................................................17 Estimation............................................................................................................................17 Quantification of Persistent vs. Transient Variation...........................................................18 Alternative models ...............................................................................................................19 DISCUSSION ...........................................................................................................................20 SUPPLEMENTARY INFORMATION ....................................................................................26 Alternative models: Estimation ...........................................................................................26 Alternative models: Results .................................................................................................28 Alternative models: Discussion ...........................................................................................29 APPENDIX ...............................................................................................................................31 LITERATURE CITED..............................................................................................................44 CHAPTER TWO: EFFECTS OF INDIVIDUAL VARIABILITY IN GROWTH ON MANAGEMENT REFERENCE POINTS: AN APPLICATION FOR LAKE TROUT (SALVELINUS NAMAYCUSH) IN LAKE SUPERIOR ................................................................51 ABSTRACT ..............................................................................................................................51 INTRODUCTION ....................................................................................................................53 METHODS ...............................................................................................................................58 Overview..............................................................................................................................58 Individual variation in growth ............................................................................................59 Simulated growth.................................................................................................................61 Calculating weights (𝛾) for each growth morph .................................................................63 Simulated selectivity and maturity schedules ......................................................................63 v Calculation of age-specific averages from known data ......................................................64 Yield-per-recruit analyses ...................................................................................................68 Spawning stock biomass-per-recruit analyses ....................................................................69 Reference points ..................................................................................................................70 RESULTS .................................................................................................................................71 Simulated growth.................................................................................................................71 Calculation of 𝑁𝑎 , 𝐶𝑎 , 𝑠𝑒𝑙𝑎 , 𝑊𝑎𝐻 , 𝑊𝑎𝑠𝑝 , and 𝑝𝑚𝑎𝑡𝑎 ..........................................................71 Standard vs. true SSBR and YPR curves .............................................................................72 Standard method estimates vs. true reference points ..........................................................73 DISCUSSION ...........................................................................................................................74 APPENDIX ...............................................................................................................................82 LITERATURE CITED..............................................................................................................95 vi LIST OF TABLES Table 1: Parameters defined and used in our von Bertalanffy model ............................................32 Table 2: Simulation cases and their assumptions of error structure ..............................................33 Table 3: AIC values for models considering different underlying random effect distributions ....33 Table 4: Parameter estimates for Model D, the best-fitting model to the data ..............................34 Table 5: Summary statistics of the frequency distribution of length-at-age for each case and each age (4, 15, and 40)..........................................................................................................................35 Table 6: Parameter estimates for the male-specific VGBM model ...............................................36 Table 7: Parameter and variable names, symbols, descriptions, and units (where appropriate) ...83 Table 8: von Bertalanffy parameter values used............................................................................85 Table 9: Biological reference points ..............................................................................................86 vii LIST OF FIGURES Figure 1: Location of the six populations of lake trout sampled from Lake Superior ...................37 Figure 2: Structure of the hierarchical von Bertalanffy growth model ..........................................38 Figure 3: Process error (transient variation) and observation error in the hierarchical model ......39 Figure 4: Biphasic growth..............................................................................................................39 Figure 5: Median estimates and 95% probability intervals for values of VBGM parameters .......40 Figure 6: Frequency distributions of length-at-age for ages 4, 15, and 40 ....................................41 Figure 7: Frequency distributions of length-at-age for ages 4, 15, and 40 using the biphasic model..............................................................................................................................................42 Figure 8: Frequency distributions of length-at-age for ages 4, 15, and 40 ....................................43 Figure 9: Individual length-at-age of Lean Lake Trout .................................................................87 Figure 10: Growth parameter values used to simulate growth morphs .........................................88 Figure 11: Individual-specific selectivity patterns and individual-specific growth patterns .........89 Figure 12: Numbers-at-age ............................................................................................................90 Figure 13: Catch-at-age..................................................................................................................90 Figure 14: Catch-at-age as a proportion of numbers-at-age (i.e. exploitation rate).......................91 Figure 15: Selectivity-at-age ..........................................................................................................91 Figure 16: Proportion mature-at-age..............................................................................................92 Figure 17: Weight-at-age (spawning & harvest) ...........................................................................92 Figure 18: Estimated YPR curves ..................................................................................................93 Figure 19: Estimated SSBR curves ................................................................................................93 Figure 20: Average weight-at-age at different levels of F.............................................................94 viii INTRODUCTION Fish growth within a population has a considerable influence on the dynamics of the system. For example, the size composition of a population – the distribution of fish sizes, sizes summarized by age, etc. – depends on how fast fish grow, how large they get, and how these parameters interact with processes like spatial and temporal variation or environmental stochasticity. Estimating fish growth rates is a complex process, typically requiring fitting statistical models to growth data and making assumptions about a growth function and error structure. Growth models often are incorporated into larger population models that provide information about stock status and condition, and thus sustainability of the stock, under exploitation by fishing. Correctly characterizing growth, then, is of major importance to both improving understanding of fishery systems and applications of modeling to fishery management. One key assumption often made when fitting growth models to length-at-age data from a fish population is that the mean length-at-age is representative of the population. More specifically, this assumes that the mean length of an age-a fish is true for all individual fish of age-a, and that estimation of average population parameters made with this assumption are true. However, individual fish vary in their growth, and ignoring individual variation in growth models has been shown to cause bias in estimation of population parameters, predictions of individual growth trajectories, and even estimation of biological reference points that inform management (DeAngelis et al., 1993; Kraak et al., 2019a; Peacor et al., 2007; Pilling et al., 2002; Sainsbury, 1980; Vincenzi et al., 2014). As an example, individual variation in size-at-age mediates how susceptible a fish is to mortality, especially when that source of mortality is size-dependent. Fishing gear is one way to regulate removal of fish from a population and is often size based: the mesh size of a gill net restricts what size of fish is caught by the net. Therefore, when a 1 population of fish is subjected to size-selective mortality, some individuals may be removed from the population at a different rate than others based on their growth. As such, these variable outcomes can have repercussions for estimates of abundance over time, projected harvest, and other population-level analyses. Estimating the magnitude and character of individual variability is one area of focus for improving growth models. Another area of interest is the source of the variation: there are two broad categories, referred to as persistent and transient [sensu (Webber & Thorson, 2016)]. Persistent variation stems from intrinsic differences in individual fish, like heritable traits, that affect how a fish grows over its entire lifespan. Transient variation, in contrast, is external and can stem from environmental conditions, year-to-year fluctuations, etc. Historically, it has been difficult to disentangle these sources of variation and in addition, measurement error can contribute to noise in growth data. If persistent variation is the dominant source, it could indicate that individual variation may be in some way heritable, potentially resulting in evolutionary consequences for a given population. If transient variation is the dominant source, it could support the importance of including environmental data in growth models. While incorporating individual variability improves the accuracy of growth models, it is also important to consider what the output of growth models are used for. One such application is a suite of models that calculate the total yield or spawning stock biomass an individual recruit (fish joining a population) will produce over its lifespan. These models, jointly termed per-recruit models, incorporate information about the average weight of fish at different ages (Quinn & Deriso, 1999b), and are used to calculate biological reference points that are used in fishery management. There have been some studies on the influence of individual variability on these models and reference points, but it is still relatively unknown how these results might be 2 generally applicable. Further, while some focus on temporal, spatial, or environmental variability in growth (Lowerre-Barbieri et al., 1998; Miller et al., 2018; Parma & Deriso, 1990; You-Gan Wang & Thomas, 1995), it is rarer to find models that explicitly model persistent and transient individual variation. Thus, for any species under exploitation that is managed using per-recruit models, additional study is required to explore whether ignoring variation introduces bias that would impact the sustainable management of the stock. Lake trout (S. namaycush) are a top piscivore in the Laurentian Great Lakes and, in Lake Superior, have a history of low population levels due to aquatic invasive species as well as supplementation from stocked fish. Several population units of lake trout in Lake Superior fall under the 1838 treaty waters management framework and thus are evaluated regularly using a statistical catch-at-age model that uses mean weight-at-age. Lake trout are known to display individual variation in growth trajectories and thus are a good application for understanding how this could influence parameter estimates and the performance of per-recruit models. In this thesis, we build on a body of research that has been characterizing the effects of ignoring individual variation in growth models. In Chapter 1, we develop a statistical model that accounts for among-individual and among-population variability for six populations of lake trout (Salvelinus namaycush) in Lake Superior. We explain the model structure, assumptions therein, and provide estimates of individual variability for the system. We also simulate from the model to assess whether the estimated variation in growth is attributable to persistent or transient variation. In Chapter 2, we take the parameter estimates from Chapter 1 and simulate a population of fish displaying individual variability in growth parameters. We then subject this population to size-selective mortality and, using parameterizations from the existing management statistical catch-at-age model, evaluate the performance of ‘standard’ per-recruit 3 models (that ignore individual variability) and ‘true’ per-recruit models (that incorporate individual variability). Our study provides quantitative measures of the magnitude and direction of individual variability for this system, show that persistent error dominates much of this observed variation, and finally, show how spawning stock and yield estimates can be biased when the variation is ignored. 4 CHAPTER ONE: MODELING INDIVIDUAL VARIABILITY IN LAKE TROUT GROWTH AUTHORS Elizabeth Stebbins1 James R. Bence1 Travis O. Brenden1 Michael J. Hansen1 1 Quantitative Fisheries Center, Michigan State University ABSTRACT Fish exhibit varying growth trajectories both within and among populations, which can interact with mortality to influence variation in size-at-age and population dynamics. Here, we developed and applied mixed-effects, hierarchical growth models with individual- and population-specific growth parameters to back-calculated length-at-age data from individual Lake Superior populations of lean morphotype lake trout (Salvelinus namaycush). Model estimates were used to simulate fish growth to determine to what extent variation in length at age was attributable to these sources. Models based on a von Bertalanffy growth function fit observed data better than those that assumed a biphasic growth function; variation in length-at-age for lean lake trout from Lake Superior was dominated by persistent sources of variation, which are intrinsic to the individual. Estimated among-individual variation exceeded estimated among-population variation in growth parameters. Implicit to our interpretation is the assumption that Lake Superior lake trout populations are not currently experiencing substantial size-selective mortality, such that estimated length-at-age patterns indeed reflect differences in growth rather than population-specific, size-selective mortality. This was deemed a reasonable assumption given the study system, but should be an important consideration in future applications of this type of modeling for other populations 5 INTRODUCTION The estimation of fish growth parameters is necessary for inferring the condition and dynamics of populations and can provide beneficial information for regulating fisheries that exploit them (Pauly, 1980). Growth, often represented as size-at-age, interacts with processes like size- selective mortality and influences or reflects population age structure, health, and maturation rates within a system of interest (Ricker, 1969a). Growth also reflects environmental conditions (e.g., prey availability) and is an important component in the development of fishery management objectives, such as size-based harvest regulations (Kristiansen & Svåsand, 1998). Many population models used by fishery managers to assess fishery exploitation levels, set management reference points, and evaluate harvest levels in relation to reference points, incorporate information about fish size at specific ages and model interactions between growth dynamics and other population-level processes (Quinn & Deriso, 1999b). Because incorrectly characterizing growth can influence the output of these broader assessments, substantial efforts have been directed at growth modeling to more appropriately model growth in different systems. For example, the biphasic growth model accounts for shifts in energy allocation at sexual maturity and appears to better describe growth in some systems (Lester et al., 2004). A common assumption of population models is that mean size-at-age is representative of an age class; typically, there is no accounting for how individual growth variability affects variability in fish size at a given age even though this variability can be substantial (Pilling et al., 2002). Although ignoring variability in growth simplifies model structure, model results may be biased. Population-level processes like mortality are influenced by variation in size-at-age and individual variability in growth (DeAngelis et al., 1993; Peacor et al., 2007). For example, individual size mediates mortality for fishes through processes like predation or fishing (Hansen et al., 1997; 6 Tsehaye et al., 2014). Because growth informs models that are used in management, the circumstances under which individual growth variability can be ignored must be clarified. Relying only on mean size-at-age may introduce bias into the calculation of biological reference points that are foundational for certain management actions (Kraak et al., 2019b; Punt et al., 2006). Growth estimates calculated from mean size-at-age have also been found to be less accurate for predicting individual growth trajectories than those calculated from models that account for individual variation (Vincenzi et al., 2014). When growth models are fit to either length-at-age data using models assuming no among-individual variation in growth, the estimated model parameters and reconstructions of growth can be biased (Sainsbury, 1980). Therefore, accounting for individual variation when characterizing or making predictions about growth may improve the accuracy of estimates and subsequently any analyses based on those estimates. Studies incorporating among-individual variation in growth are limited worldwide, including in North America’s Laurentian Great Lakes region. Variation in size-at-age can be broadly categorized as transient or persistent. Persistent variation comes from individual-specific characteristics that influence how a fish grows over its lifespan, possibly stemming from heritable attributes or early-life environmental conditions that manifest over the course of the individual’s life. Transient variation stems from year- and individual- specific differences in growth of a fish, potentially capturing environmental variation or year effects (Morrongiello et al., 2015; Pfister & Stevens, 2002; Webber & Thorson, 2016) . An important characteristic of transient variation is that, because it is individual- and year-specific, it is not a direct corollary of year effects like those modeled in Weisberg et al. (Webber and Thorson, 2016; Weisberg et al., 2010). Systems where among-individual growth variability is dominated by persistent variation may be more susceptible to processes like size-selective mortality, such as 7 a population that becomes dominated by slow-growing individuals due to selective harvest of faster-growing individuals (Biro & Post, 2008; Webber & Thorson, 2016). Quantifying how persistent and transient components of growth variation contribute to observed variation in size- at-age can be challenging, but mixed-effects modeling offers a potential solution (Webber and Thorson 2016). For Antarctic toothfish (Dissostichus mawsoni), nearly half of the variation in length-at-age was attributable to transient variation (Webber and Thorson, 2016), with the rest allocated to observation error (persistent and sex-specific variation were not detected; Webber and Thorson, 2016). However, just because transient growth was the primary source of variation for this particular species does not mean persistent variation is never a significant source of variation in fish size at age. In systems where persistent variation is dominant and is tied to heritability, evolutionary trajectories of the population could be affected (Biro & Post, 2008; Wolf & Weissing, 2012); additionally, the magnitude of persistent variation in growth could determine the extent to which size-selective mortality affects the size structure of the stock, and subsequently introduce bias into management reference points based on growth parameters (Parma & Deriso, 1990). A key point here is that the characteristics of among individual variation in growth has importance in its own right, not merely as nuisance parameters that allow for better estimation of mean growth patterns. While there are increasing applications of hierarchical growth models for fish, including lake trout that allow for among individual variation in growth (Chavarie et al., 2017, 2019; Hansen, Nate, Chavarie, et al., 2016), many of these do not actually report on the variance parameters, meaning the published literature provides limited understanding of the nature of growth variation. North America’s Laurentian Great Lakes are home to multiple fish species that are commercially and recreationally important. Most of the population assessment models used to manage Great 8 Lakes fisheries, like catch-at-age models, assume mean length-at-age is representative of the assessed population. Similarly, models that have calculated prey fish consumption in the Great Lakes have assumed that all fish at a given age were the same size (Tsehaye et al., 2014; He et al., 2015). The history of the Great Lakes ecosystem has been characterized by multiple invasions of aquatic nuisance species [e.g., sea lamprey (Petromyzon marinus); alewife (Alosa pseudoharengus)] that have drastically altered fish community composition. Historically, lake trout (Salvelinus namaycush) were the predominant piscivore throughout the Great Lakes, except Lake Erie, but their abundances were severely reduced due to sea lamprey predation, overfishing, declining water quality, and spawning habitat destruction (Hansen, 1999). Beginning in the mid- 1900s, stocking was initiated by Great Lakes fishery agencies to rehabilitate or restore lake trout populations in each lake. A wide number of lake trout strains were stocked to identify the strains most likely to restore self-sustaining populations. To date, Lake Superior is the only Great Lake where lake trout rehabilitation efforts have successfully led to self-sustaining populations to the point where stocking has been deemed no longer needed (Muir et al., 2012). Even though stocking is no longer conducted in Lake Superior, the genetic makeup of wild lake trout may reflect this complex stocking history (“The Lake Charr Salvelinus Namaycush: Biology, Ecology, Distribution, and Management,” 2021). Lake trout in Lake Superior thus represent a case study of a population with a history of disturbance, low abundance, and rehabilitation efforts involving stocking fish with varied genetic composition. Measuring individual variability and understanding its sources could be an important addition to characterizing lake trout populations in Lake Superior. The aim of this study was to quantify the magnitude of among-individual variability in growth as well as the relative contribution of transient and persistent variation to lake trout in Lake 9 Superior, with potential applications to other important fish species. We used a mixed-effects, hierarchical model to estimate transient and persistent growth, among individuals and populations, for a length-at-age dataset of six populations of lake trout in Lake Superior. We also estimated correlation among growth parameters among individuals and populations. We build on a body of work that has used mixed-effects models to partition and estimate sources of variability in fish growth (Hart & Chute, 2009; Ogle et al., 2012; Vigliola & Meekan, 2009; Webber & Thorson, 2016; Weisberg et al., 2010). We hypothesized that a multivariate normal distribution for von Bertalanffy parameters (with necessarily positive parameters on the log-scale) would account for greater variability than univariate (non-correlated) distributions, and that L∞ and K would be negatively correlated among populations and individuals. Further, we hypothesized that incorporating separate persistent and transient variation would allow us to understand their relative contribution to length-at-age for lake trout. 10 METHODS Data Collection Lake trout were collected from six populations in Lake Superior in 2002, 2003, 2004, 2006, 2007, 2013, and 2014 (Figure 1). Sex, length, and weight assignment processes, as well as otolith processing, are as described in Hansen et al. (2016b). Four morphotypes of lake trout are generally recognized as existing in Lake Superior, but visual discrimination of these morphotypes can be unreliable, with the lean morphotype being the most reliably identified visually (Muir et al., 2014). To reduce uncertainty, only fish identified as lean morphotype trout were used in this analysis (n=410). After otoliths were cleaned, increment width for otolith annuli were measured. Length at age of each individual fish was back-calculated using the biological intercept model that assumes a linear relationship between otolith radius and fish length with a biologically determined intercept (Campana, 1990). Biological intercept values were obtained by fitting the biological intercept model to length- and otolith radius-at capture (Vigliola & Meekan, 2009). The von Bertalanffy Growth Model for an Individual Fish The aim of the hierarchical growth model was to estimate growth function parameters for each individual fish i belonging to population p. The von Bertalanffy growth model (VBGM) is the most commonly used growth model in fisheries science (Flinn & Midway, 2021). The Fabens (1965) parameterization of the VBGM was used to model growth increments as a function of L∞ = asymptotic length, K = Brody growth coefficient, and L1 = length at age 1: 𝐿(1) = 𝐿1 (1.1) 𝐿 (𝑎) = 𝐿 (𝑎 − 1) + ΔL(𝑎 − 1) for 𝑎 ≥ 2 (1.2) ΔL 𝑎) = (𝐿∞ − 𝐿 (𝑎) ∗ (1 − exp(−𝐾 )) for 𝑎 ≥ 2 ( (1.3) 11 That is, length at the start of the annual growth season for the first age L(1) is set equal to the estimated parameter L1 (Equation 1.1), with length at the start of the year for later ages increasing by an growth increment ΔL from the length at the start of the previous years (Equation 1.2). The ΔLs depend upon length at the start of the year, L∞, and K (Equation 1.3). Variation among individuals in VBGM parameters (L∞, K, L1) constitutes persistent among-individual variation, and the hierarchical way we modeled this variation is described in the Hierarchical Model Structure Section. Equation 1.3 was modified to allow for transient variation: ΔL(𝑎) = (𝐿∞ − 𝐿(𝑎)(1 − exp(−𝐾)) ∗ exp(𝛿(𝑎)) (2) where 𝛿(𝑎) are log-scale age-specific process errors drawn separately for each year of growth (and separately for each individual) independently from a common normal distribution: −𝜎𝛿2 2 (3) 𝛿(𝑎)~𝑁( , 𝜎𝛿 ) 2 −𝜎𝛿2 The expectation for the normally distributed log-scale process errors term is rather than zero 2 so that the expectation of the multiplicative process error exp(𝛿(𝑎)) is 1.0. Formally, the 𝛿 (𝑎) (for each growth year for each fish) are random effects. An issue with using equations 2 and 3 as a replacement for equation 1.3 is that an individual could exceed its individual asymptotic size following a year with a large positive 𝛿(𝑎), and thus need to shrink in the next year. A simple fix to this issue would be to set the growth increment to (𝐿∞ − 𝐿(𝑎)) whenever this occurred. This, however, causes estimation issues by making the model non-differentiable. We implemented a differentiable approximation to this approach by setting the true growth increment equal to the weighted average of the right-hand side of Equation 2 and (L∞-L(a)), where the weights for the former and latter rapidly shift from 1 to 0 12 and from 0 to 1 respectively as the right-hand side of Equation 2 approaches and exceeds (L∞- L(a)): ΔL(𝑎) = (𝑤2 ∗ Δ ̂𝐿(𝑎)) ∗ (𝑤1 ∗ (𝐿∞ − 𝐿(𝑎)) (4.1) 1 𝑤2 = 1 − exp (−𝛼 ∗ (0.99 ∗ 𝐿∞ − (𝐿(𝑎) + (Δ ̂𝐿(𝑎))) (4.2) w1 = 1 − w2 (4.3) ̂ (𝑎) = (𝐿∞ − 𝐿(𝑎))(1 − exp(−𝐾 )) ∗ exp (𝛿 (𝑎)) ΔL (4.4) We set 𝛼 equal to 2, selected to ensure the transition from one weight to the other was sufficiently steep to happen within a year. Hierarchical model structure Among-individual variability in the VBGM parameters was modeled using a hierarchical structure, so that in the second level mean values of log-scale parameters for each population {log 𝐿∞𝑝 , log 𝐾𝑝 , log 𝐿1𝑝 } varied around population-averaged means {log 𝐿∞∙∙ , log 𝐾∙∙ , log 𝐿1∙∙ } (Eq. 5.1), and in the first level the VBGM parameters for an individual {log 𝐿∞𝑝,𝑖, , log 𝐾𝑝,𝑖 , log 𝐿1𝑝,𝑖 } varied around their population means (Eq. 5.2, Figure 2). Variation in both levels of the hierarchical model was assumed to follow a multivariate normal distribution, with level-specific variance-covariance matrices (Equations 5.1 & 5.2).: {log 𝐿∞𝑝 , log 𝐾𝑝 , log 𝐿1𝑝 } ~ 𝑀𝑉𝑁({log 𝐿∞.. , log 𝐾.. , log 𝐿1.. }, Σ.. ) (5.1) {log 𝐿∞𝑝,𝑖 , log 𝐾𝑝,𝑖 , log 𝐿1𝑝,𝑖 } ~ 𝑀𝑉𝑁 ({log 𝐿∞𝑝∙ , log 𝐾𝑝∙ , log 𝐿1𝑝∙ } , Σ𝑝 ) (5.2) Thus, persistent variation among all individuals was a consequence of variation in mean parameter values among populations and variation around these means within populations, because these influenced a fish’s growth trajectory over the course of its life. Both population 13 means and values for each individual within a population of the VBGM were both random effects. Distribution of observed data The previous sections (The von Bertalanffy Growth Model for an individual fish and Hierarchical model structure) described how a fish was assumed to grow given its individual growth parameters and how those growth parameters varied among individuals and populations. Conditioned on the individual specific VBGM parameters and process errors, Equation 1 produced what is assumed in the model to be the true sequence of lengths at age for each fish up to the time it was collected. Observed lengths at age were assumed to have expectations equal to true lengths at age that followed a lognormal distribution: 𝜎𝜀2 2 log 𝐿̃ (𝑎, 𝑝, 𝑖 ) ~𝑁(log(𝐿(𝑎, 𝑝, 𝑖 )) − ,𝜎 ) 2 𝜀 (6) Equation 6 represents the incorporation of random observation error with a bias correction term 𝜎𝜀2 of − in the mean and a standard deviation of 𝜎𝜀2 (Figure 3). 2 Estimation Template Model Builder (TMB) was used for model fitting because it is designed for highly parameterized non-linear models with random effects. TMB uses C++ code and an R software interface to compile and implement automatic differentiation, integrating efficiently over random effects via Laplace approximation (Kristensen et al., 2016). A combination of fixed and random effects were estimated using maximum likelihood estimation (Table 1). TMB provides options for efficiently parameterizing the variance-covariance matrices ∑𝑝 and ∑∙∙ in the form of standard deviations for each parameter and a set of parameters that determines the correlations among them (Kristensen et al., 2022). 14 In preliminary efforts to fit the model, models that allowed all three random effects for the parameters to be correlated at both the individual and population level failed to converge. Consequently, parameterizations for which only two of the VBGM parameters were correlated, with the third being uncorrelated with the others within a level, or for all three parameters to be uncorrelated. Thus, combinations of a multivariate normal (all three VGBM parameters correlated), bivariate normal (just L∞ and K correlated), and univariate normal (all three VBGM parameters uncorrelated) were modeled. A bivariate normal with L∞ and K did not consider L1 in any bivariate normal case. L∞ and K represent a relationship between body size and speed of growth whereas L1 only represents the first year of growth. Several analyses have found L∞ and K to have a significant relationship, both positive and negative, which could either suggest that size ranks (i.e., larger fish grow faster, smaller fish grow slower) are maintained or not maintained (Vincenzi et al., 2014). The best-fitting model based on AIC was used for subsequent analyses. Quantification of Persistent vs. Transient Variation Estimates of variation and covariation provided information on persistent and transient variation. However, different variances contribute in different ways to total variation in observed lengths given age that is not readily evaluated by simply examining estimated parameters. Consequently, to translate parameter estimates into estimates of the relative contribution of persistent and transient variation (and observation error) into observed variation among fish in length at age, datasets based on fitted models were simulated for versions of models with some sources of variation removed. Three cases were simulated: (1) a base case, with persistent and transient variation, (2) a case without transient variation, and (3) a case without persistent variation (Table 2). For each case, growth of 1,000 fish was simulated from ages 1 to 43 years 15 and the frequency distribution of length-at-age in simulated datasets was examined for selected ages. By comparing the distributions of cases two and three, which represent systems with no transient error and no persistent error (respectively), to cases one and two, which represent systems with all sources of random error and no sources of random error except observation error (respectively), we determined how each type of error influenced variation in observed length-at-age. Alternative Models: Biphasic and Sex-Specific Criticisms of the VBGM include the parameters not being biologically interpretable and that it does not account for shifts in energy allocation when fish reach maturity, which may affect growth patterns (Lester et al., 2004). Lake trout are late maturing and long-lived, and some individuals in our sample potentially displayed biphasic growth, where immature growth was linear and growth after maturity followed a VBGM pattern (Figure 4). Therefore, models based on a Lester biphasic growth function (Lester et al., 2004) were evaluated, as well as versions of the VBGM applied to male and female lake trout data separately to test for sex-specific patterns (Supplementary materials). 16 RESULTS Estimation The model that best fit the data, Model D, assumed that individual log-scale growth parameters L∞, K, and L1 were distributed as multivariate normal, population-level means for log-scale L∞ and K were defined by a bivariate normal distribution, and log-scale L1 was distributed normally and was uncorrelated with L∞ and K. This model outperformed the next best model by 8 AIC units (Table 3). Because model D fit substantially better than the alternative models, detailed results were described only for this model, and it was used as the basis for quantification of persistent and transient variation. Persistent variation was greater for all three growth parameters at the individual level than at the population level (Figure 5). Growth model parameter estimates varied considerably among individuals after accounting for population-level variability (Table 4). Back-transformed (median) estimates [95% CIs] for population-averages for the VBGM parameters were 728 mm [633, 837] for asymptotic size, L∞, 0.10 y-1 [0.08, 0.11] for the Brody growth coefficient, K, related to how fast asymptotic size was approached, and 106 mm [99, 112] for length at age 1, L1. Standard deviations of log-scale K were 0.13 [0.05, 0.36] at the population level and 0.33 [0.29, 0.38] at the individual level. Standard deviations of log-scale L∞ were 0.15 [0.08, 0.28] at the population level and 0.24 [0.22, 0.26] at the individual level. Standard deviations of log-scale L1 were 0.06 [0.03, 0.13] at the population and 0.20 [0.18, 0.22] at the individual level. The standard deviation of observation error was 0.0004 [0.0003, 0.0007]. Finally, the estimated standard deviation due to process errors was 0.2346 [0.2254, 0.2441] (Table 4). For all parameters, estimated among-individual standard deviations were larger than estimated among-population standard deviations, but K differed the most between levels and L1 differed the least. At the individual level, the estimated standard 17 deviation for log-scale K was larger than that for log-scale L∞. At the individual level, log-scale L∞ and K were negative correlated (-0.78), log-scale L∞ and L1 were positively correlated (0.02), and log-scale K and L1 were positively correlated (0.21). At the population level, log scale means for L∞ and K were positively correlated (0.63). Correlations were opposite for population-level and individual-level relationships of L∞ and K. At the population level, L∞ and K were positively correlated (0.630) and at the individual level, negatively correlated (-0.78). Correlations between L1 and other VBGM parameters were not estimated at the population level, but at the individual level, L1 and L∞ were positively correlated (0.023) and L1 and K were positively correlated (0.213). Quantification of Persistent vs. Transient Variation Simulated distributions of length-at-ages 4, 15, and 40 were similar for the model with both persistent and transient errors (case one) and for the model with only persistent error (case two) (Figure 6). Density distributions for the model with only transient error (case three) were, for all ages, narrower than for models with both persistent and transient errors (case one) and for only persistent error (case two) but were very narrowly distributed around a mean of 714.9 mm (Figure 6, Table 5). Standard deviations of length-at-age for the case with both persistent and transient error, as well as just persistent (cases one and two) increased with age, but for the case with just transient error (case three) decreased. Going from age 15 to age 40, the variance in length-at-age for only transient error (case three) dropped from 16.59 to 3.46. Mean of age-specific distributions were relatively constant among cases, but not among ages (Table 5). 18 Alternative models Persistent and transient variation was difficult to estimate for the biphasic model and sex-specific VBGM growth function (see Supplementary materials), and did not suggest markedly different conclusions than were reached with the base VBGM. 19 DISCUSSION An implicit assumption we made in our model was that variation in observed size-at-age was not due to size-selective mortality, which can interact with growth (Ricker, 1969a). Based on levels of mortality for lake trout in Lake Superior, we assumed that size-selective mortality was not sufficient over time to have substantially altered distributions of growth parameters of survivors (e.g., so fish with a relatively high L∞ or K were rarer at certain ages) (Lenart & Caroffino, 2019). We found that among-individual variability in the VBGM parameters for a population of lake trout in Lake Superior exceeded the variability in population means for each parameter. These results suggested that lake trout, at least in Lake Superior, exhibited significant among- individual variation in growth among fish collected from localized populations. Lake trout morphotypes exhibit variation in growth parameters across their range (Chavarie et al., 2017; Hansen et al., 2016a, 2016b, 2012). Some physiological traits that differ among lake trout morphotypes are heritable (Wellband et al., 2021), and morphotype to some extent determines growth parameter estimates. Our study focused on lean morphotype lake trout with location- based populations but similarly shows a population-level effect on growth parameter estimates. This reflects the complexity of discerning the source of individual variation in growth parameters, especially for species like S. namaycush, where resource use as well as genetic differences could be contributing to persistent variation. Our model could easily be applied to datasets with additional populations from more than one lake. This would be particularly illuminating if done with lake trout to further elucidate whether our results were specific to Lake Superior or if variation within localized populations exceeds that of variation among populations across the species’ range. 20 Our study partitioned variation in growth into persistent and transient sources; however, disentangling intrinsic and extrinsic variation is not a novel concept – for example, recent approaches partition individual variation into similar categories (Stawitz et al., 2015, 2019; Thorson & Minte-Vera, 2016; Vincenzi et al., 2014; Weisberg et al., 2010). Our model treats transient (i.e., extrinsic) error as individual- and year-specific to explicitly model the effects of extrinsic sources on individual-specific variation from year-to-year. There are two key reasons to focus on understanding the magnitude and nature of persistent variation: if measured intrinsic variation were heritable it could have significant impacts on the evolutionary trajectory of a population and represent genetic diversity in a population (Biro & Post, 2008; Wolf & Weissing, 2012). In this case, estimating the magnitude could potentially be connected to the magnitude of evolutionary effects. Further, systems under exploitation where persistent variation is dominant could exhibit bias in length-based stock assessments and shifts in population structure (Parma & Deriso, 1990). Estimating the magnitude of persistent variation, and its relative contribution to individual variation in growth, across a range of populations, species, and locations, would further elucidate how persistent variation ties into these broader dynamics and help direct further inquiry. With the development of a range of alternative models that allow for among individual variation as part of hierarchical analyses, models that account for among individual variation are being increasingly applied. Unfortunately, often this variation is treated as being nuisance parameters so that in the majority of such applications the variance parameters are not even reported. We urge the regular reporting of these variance parameters so that we can build up the information needed to better reach generalizations regarding growth variation. To understand the consequences of ignoring persistent (among-individual and among- population) variation in growth, additional analysis is required, a topic we return to in Chapter 2. 21 Variability in fish size can interact with management analyses such as spawning stock biomass, wherein fisheries are managed to conserve some level of spawning stock in a population (Quinn & Deriso, 1999b). The number of eggs deposited by a sexually mature female lake trout depends on the size of the female (Hansen et al., 2021), and individual large, fecund female fish play an important role in fish population dynamics (Hixon et al., 2014). Spawning stock biomass per- recruit calculations often used to set fishing mortality reference points generally are calculated using mean size-at-age and assume that size selective fishing mortality will not influence these means (Caroffino & Lenart, 2000; Lenart & Caroffino, 2019). We demonstrated individual variability in growth of lake trout and showed that this variation was largely due to intrinsic differences in VBGM parameters, which could interact significantly with predictions of spawning stock biomass per recruit or similar mean weight-at-age-based analyses. For example, applying size-selective mortality specific to individual fish with different growth parameters and associated length-at-age schedules would allow us to understand how this influences population structure and potential evolutionary consequences. Allowing individuals to vary in their fishing mortality risk can significantly influence subsequent yield-per-recruit estimates, thereby suggesting that understanding individual-specific fishing vulnerability could be a productive line of inquiry for our model (Hart, 2001). Individual variability, modeled through individuals being assigned different growth parameters, interacted significantly with population dynamics through stock-recruitment functions, heritability, and different fishing regimes, thereby altering the distribution of growth parameters in later generations (Martínez- Garmendia, 1998). Further, the potential risks of mischaracterizing age and size structure of populations due to incorrect simplifications are high for management (Punt et al., 2017; Webber & Thorson, 2016). 22 Log-scale variation in growth increments due to transient error had a standard deviation of 0.234, which suggests that external sources of variation from year-to-year influence how an individual grows. This magnitude cannot be directly compared to variation in the VBGM parameters that lead to persistent variation because they are on different time scales, and variation in VBGM parameters passes through the nonlinear growth model to produce variation in increments. Our results contrast with another study that found transient error was important whereas persistent variation played little role in determining variation in length-at-age of Antarctic toothfish (Dissostichus mawsoni) (Webber & Thorson, 2016). This difference suggests that the role of different types of variation may depend upon specific characteristics of fish species and/or the surrounding environment, although too few examples exist to generalize. The nature of the study system may be important for the significance of our findings. Lake Superior is deeper and larger and on average shows lower interannual variability in surface temperature than the other Great Lakes. Similarly, temporal variability of the environment certainly differs among inland lakes. Applying our model to populations of lake trout in other lakes might improve understanding as to patterns of variability that uncover generalities and explanations for system specific results. The individual variation in L1, which we categorized as representing persistent variation, could also be interpreted as a combination of persistent and transient variation. The length at age 1 and variation thereof represents growth during just the first year, which could in some cases reflect conditions like density dependence, temperature, or other year-specific fluctuations. For example, density during the first year of life and year-of-birth effects were significant predictors of L∞ and K for two populations of Marble trout (Salmo marmoratus) (Vincenzi et al., 2014). Theoretically, density dependence or environmental conditions in the first year of life may have influenced our estimates of L1 variation, which we treated as persistent. In our model, no method 23 can separate L1 variation associated with growth over the lifespan (persistent variation) from what is associated with early life conditions. However, we saw no year-specific patterns in our estimates of transient variation across the dataset. The influence of initial length on persistent and transient growth variation may be too difficult to estimate accurately from common fisheries data like catch time series (Parma & Deriso, 1990). We believe our model can approximate the relative contribution of persistent and transient variability, although we acknowledge that lack of access to similar high quality data in other systems would limit the application of our method. We encountered convergence issues when we attempted to treat all three VBGM parameters as correlated on a multivariate lognormal distribution at both individual and population levels. We ultimately tested several different assumptions of correlation among growth parameters, both among individuals and among populations, to find the best-fitting model that converged. We encountered this issue with both the VBGM and biphasic model. The best fitting VBGM model treated L∞, K, and L1 as correlated among individuals, and only L∞ and K as correlated among populations. This suggests that data from this system did not support an assumption of correlation among all parameters at the population level. Past studies allowed one parameter to vary (James, 1991; Kimura et al., 1993; Kirkwood & Somers, 1984; Punt et al., 2006; Wang, 1998), while other studies allowed L∞ and K to covary (Hart & Chute, 2009; Pilling et al., 2002; Xiao, 1994; Zhang et al., 2009), but this correlation has been both positive and negative depending on the study system (Eveson et al., 2007; Ortiz de Zarate & Babcock, 2015). A positive correlation between L∞ and K may indicate that growing large early in life leads to larger asymptotic size, whereas a negative correlation may indicate an energetic tradeoff that reduces asymptotic size of individuals that grow large at an early age (Vincenzi et al., 2014). Our results suggest that covariation of L1 with L∞ and K was negligible at the population level, 24 although the relatively small number of populations in our study (six) and the associated lack of statistical power may have contributed to difficulties estimating among-population variation or covariation among VBGM parameters. At the individual level, however, we found that the correlation between L1 and K was an order of magnitude stronger than that of L1 and L∞. Our results suggest a positive relationship between L1 and K, but we are not aware of other studies that used L1 or allowed L1 to vary among individuals. However, further studies using a trivariate normal for VBGM parameters could shed light on whether this relationship is common among species or populations of lake trout. We found persistent individual- and population-level variability in growth of lake trout in Lake Superior and provided estimates of that variation and its sources. Our hierarchical model, based on a von Bertalanffy growth function, performed better than alternative models like the biphasic, which has been shown to work for lake trout. We were able to estimate measures of variability for all three von Bertalanffy parameters and their covariation and note that this variability is significant. We acknowledge that individual variability in growth, when ignored, can bias population-averaged parameter estimates and that this bias can be propagated into management reference points, although many models do not account for it. Our study develops a complex, hierarchical model that can be fit to data from multiple populations and presents measurements of spread for growth variation that can be helpful to those concerned with the Laurentian Great Lakes, but also to any system where individual variability is suspected. 25 SUPPLEMENTARY INFORMATION Alternative Models: Estimation We adapted a biphasic growth model that predicts growth as linear during the immature stage until an inflection point (intended to be age of maturity, represented by T), when growth switches to a VBGM curve. The model was developed by Quince et al. (2008a,b) and later promoted by Lester et al. (Lester et al., 2014; Quince et al., 2008a, 2008b). Our model grows an individual fish i at age a, beginning with length at age 1: 𝐿(1, 𝑖 ) = 𝐿1 𝑖 (S1) Growth increments were calculated and added at each time step, like the base model (Equation S2). 𝐿(𝑎, 𝑖 ) = 𝐿(𝑎 − 1, 𝑖 ) + ΔL(𝑎, 𝑖 ), 𝑎 > 1 (S2) While fish were immature (a < T), the growth increment was dominated by linear growth and defined by hi, the average annual immature growth rate (S3). While the fish was mature (a ≥ T) the growth increment was dominated by VBGM growth (S4). Δ̇𝐿(𝑎, 𝑖 ) = ℎ𝑖 (S3) Δ̈𝐿(𝑎, 𝑖 ) = (𝐿∞𝑖 − 𝐿(𝑎 − 1, 𝑖 )(1 − exp(−𝐾𝑖 )) (S4) To avoid the likelihood being a non-differentiable function of the inflection point T, growth changes gradually from linear to VBGM over a one-year period. This is achieved by modeling growth increments (S5) as a weighted sum of the linear and VBGM increments, with the weights (Equation S6 and Equation S7) summing to 1 and the weight for the linear increment going from near 1 to near 0 during the year of maturation. As a approaches T, the linear weight decreases according to a logistic function (Equation S6). Δ𝐿(𝑎, 𝑖 ) = (Δ̇𝐿(𝑎, 𝑖 ) ∗ 𝑤1 + Δ̈𝐿(𝑎, 𝑖 ) ∗ 𝑤2 ) ∗ exp (𝛿(𝑎, 𝑖 )) (S5) 26 1 (S6) 𝑤1 = 1 + 𝑒𝑥𝑝(−3(𝑎𝑖 − 𝑇𝑖 )) 𝑤2 = 1 − 𝑤1 (S7) The individual-specific parameters we modeled for the biphasic growth model were g, h, and L1, representing the ratio of energetic investment in gonads to somatic mass, the annual immature growth rate, and length at age 1 respectively (Lester et al., 2014). The individual-specific parameters gi and hi can be transformed into VBGM parameters 𝐿∞𝑖 (Equation S8) and 𝐾𝑖 (Equation S9). 3ℎ𝑖 (S8) 𝐿∞𝑖 = 𝑔𝑖 𝑔𝑖 (S9) 𝐾𝑖 = ln (1 + ) 3 We modeled individual-specific values for T as a linear function of hi (Equation S10), incorporating a slope parameter as β. 𝑇𝑖 = (𝛽 ∗ ℎ𝑖 ) (S10) We modeled g, h, and L1 on a lognormal scale using the same hierarchical structure as for the base model (Figure 2). Assumptions about the distribution of process errors and observation errors were the same as the base model, and similarly to the base model we assume that the population means and individual specific growth model parameters on a log-scale followed multivariate normal distributions: {log 𝑔𝑖 , log ℎ𝑖 , log 𝐿1 𝑖 }~𝑀𝑉𝑁({log 𝑔𝑝 , log ℎ𝑝 , log 𝐿1 𝑝 } , Σ𝑝 ) (S11) {log 𝑔𝑝 , log ℎ𝑝 , log 𝐿1𝑝 } ~𝑀𝑉𝑁({log 𝑔.. , log ℎ.. , log 𝐿1.. }, Σ.. ) (S12) 27 Alternative Models: Results The biphasic model produced estimates for all parameters when only among-individual variation was included in the model. When the second level of the hierarchical model, which models among-population variation for the three growth parameters, was included, the model failed to converge regardless of correlation structure for the growth parameters. This suggests that the data lack power to support among-population variation in the biphasic growth parameters. The biphasic growth model, when only among-individual variation between log-scale parameters was included, produced an AIC 1400 larger than the best-fitting VBGM model. We fit the VBGM to female and male lake trout separately. The model failed to converge when run with only female lake trout data, and this is likely not due to sample size (nF = 203, nM = 205). When we subsampled the original dataset to include only males or only females, this also could have selectively pulled from certain populations, collection ages, etc., resulting in a dataset that had some populations represented significantly differently than in the combined-sex dataset. The male-specific VBGM model parameter estimates were similar to those of the combined-sex parameter estimates (Table 6). Specifically, the variances of the individual and population-level parameters only differed slightly. The estimated mean of the log L∞ parameter was smaller for male-specific model (706.27mm compared to 727.78mm), being the only major difference. The standard errors also showed negligible differences from the combined-sex model. Examination of both simulations from both biphasic model and the male-specific VBGM model showed that persistent variation contributed more to variation in length-at-age than did transient variation, supporting our results from the best-fitting VBGM model (Figures 7 and 8). 28 Alternative models: Discussion Our biphasic model failed to converge when among-population variation was included in the model. Quince et al. applied their biphasic model to lake trout data and succeeded in fitting it (Quince et al., 2008b). There are a few key differences between our model and the biphasic model developed by Quince et al. Ours is individual-based and allows for variation at the population level, which even for the VBGM caused issues with convergence, possibly out of overparameterization. We also calculated T, or age at maturity, as a linear function of h (immature growth rate). Long-lived, late-maturing species are good candidates for the biphasic model and particularly in the case of lake trout, where the biphasic model has been shown to outperform the VBGM (Quince et al., 2008b), it was necessary to consider. However, in this case the biphasic model was unable to support both among-individual and among-population variation where the VBGM was. The VBGM has additional benefits because of its long history of use across species and taxa, which allows us to compare parameter estimates across studies (Flinn & Midway, 2021). AIC values showed that a VBGM model with two levels fit better than a biphasic model with only among-individual variation. Further, the male-specific model did not estimate significantly different parameters compared to the combined-sex model. This suggests that the VBGM model is a better fit to the data. Certainly, in this instance the VBGM model is better suited to estimating parameters in a hierarchical framework to longitudinal data. To our knowledge there are few studies attempting to fit a biphasic model to individual growth data, and further attempts to do so would provide better information on whether the VBGM model better describes individual lake trout growth, or how much the quantity and structure of the data influence the comparison. 29 Both the biphasic and male-specific VBGM models showed persistent variation contributing more to variation in length-at-age in simulated datasets. This supports our earlier results and shows that they are not specific to the best fitting VBGM model, but also apply to alternative growth functions and a subset of the data. 30 APPENDIX 31 APPENDIX Table 1 Parameter Symbol Description Observation Error 𝜎𝜀 Observation error for length increments SD 𝜎𝐿∞ Standard deviation of individual-level 𝑖 multivariate lognormal distribution L∞ SD (individual) component for L∞ 𝜎𝐾𝑖 Standard deviation of individual-level multivariate lognormal distribution K SD (individual) component for K 𝜎𝐿1 Standard deviation of individual-level 𝑖 multivariate lognormal distribution L1 SD (individual) component for L1 𝜎𝐿∞𝑝 Standard deviation of population-level bivariate lognormal distribution L∞ SD (population) component for L∞ 𝜎𝐾𝑝 Standard deviation of population-level bivariate lognormal distribution K SD (population) component for K 𝜎𝐿1𝑝 Standard deviation of population-level L1 SD (population) lognormal distribution component for L1 𝐿∞.. Metapopulation (Lake Superior) mean for log L∞ mean L∞, the von Bertalanffy parameter for (metapopulation) asymptotic length log K mean 𝐾.. Metapopulation (Lake Superior) mean for (metapopulation) K, the Brody growth coefficient 𝐿1.. Metapopulation (Lake Superior) mean for log L1 mean L1, the von Bertalanffy parameter for (metapopulation) length at age 1 𝜎𝛿 Standard deviation for year- and age- specific process error, here representing Process Error transient error Table 1: Parameters defined and used in our von Bertalanffy model. Thetas 1-4 are the parameter used by TMB to create the correlation matrices for the multivariate and bivariate distributions: Thetas 1-3 are used in creating the matrix at the individual level and Theta 4 is used in the matrix at the population level. Estimates for the theta values are not individually interpretable and are thus not reported, but the correlations between parameters they calculate are reported in Table 4. 32 Table 1 (cont’d) Parameter Symbol Description Ѳ1 Parameter used in determining correlation Theta 1 among individual log L∞ and log K means Ѳ2 Parameter used in determining correlation Theta 2 among individual log L∞ and log L1 means Ѳ3 Parameter used in determining correlation Theta 3 among individual log L1 and log K means Ѳ4 Parameter used in determining correlation Theta 4 among population log L∞ and log K means Table 2 Case Transient error Persistent error (𝝈𝟐𝜹 ) ∑ 𝑝 and ∑. . 1 YES YES 2 NO YES 3 YES NO Table 2: Simulation cases and their assumptions of error structure. Case 1, our base case, and Case 4, serve as comparison for Cases 2 and 3, which represent regimes with or without persistent and transient variation. Table 3 Model Individual-level distribution Population-level distribution AIC A BVN BVN -73133.2 B UVN UVN C MVN MVN D MVN BVN -73180.5 E BVN MVN F BVN UVN -73124.8 G MVN UVN -73172.5 H UVN MVN I UVN BVN Table 3: AIC values for models considering different underlying random effect distributions. AIC cells with a diagonal line and shading indicate the model did not converge. All models use the Fabens parametrization of the von Bertalanffy growth model. Here BVN = bivariate lognormal with L∞ and K, covarying and L1 uncorrelated; MVN = multivariate lognormal with L∞, K, and L1 covarying; UVN = univariate lognormal with L∞, K, and L1 all uncorrelated. 33 Table 4 Symbol Estimate Estimate Parameter (log) (numerical) SE (log) CI LB CI UB Observation Error 𝜎𝜀 -7.74 0.0004 0.24 0.0003 0.0007 L∞ deviations 𝜎𝐿∞ 𝑖 (individual) -1.43 0.24 0.05 0.22 0.26 K deviations 𝜎𝐾𝑖 (individual) -1.1 0.33 0.07 0.29 0.38 L1 deviations 𝜎𝐿1 𝑖 (individual) -1.61 0.2 0.04 0.18 0.22 L∞ deviations 𝜎𝐿∞𝑝 (population) -1.93 0.15 0.32 0.08 0.28 K deviations 𝜎𝐾𝑝 (population) -2.05 0.13 0.52 0.05 0.36 L1 deviations 𝜎𝐿1𝑝 (population) -2.8 0.06 0.38 0.03 0.13 L∞ mean 𝐿∞.. (metapopulation) 6.59 727.78 0.07 632.7 837.15 K mean 𝐾.. (metapopulation) -2.34 0.1 0.07 0.08 0.11 L1 mean 𝐿1.. (metapopulation) 4.66 105.64 0.03 99.48 112.17 Process Error 𝜎𝛿 -1.45 0.23 0.01 0.23 0.24 L∞ / K correlation (individual) -0.778 0.118 L∞ / L1 correlation (individual) 0.023 0.064 K / L1 correlation (individual) 0.214 0.398 L∞ / K correlation (population) 0.630 0.812 Table 4: Parameter estimates for Model D, the best-fitting model to the data. All estimates are presented as standard deviation of the underlying distribution except the metapopulation means. Estimates are presented on the log scale they were estimated and back-transformed to the numerical scale. Confidence intervals were calculated as (ESTIMATE ± 2*SE ) for parameters estimated on the log scale and then the lower and upper bounds were back-transformed (by exp). Correlations were not estimated on a log scale but instead determined from estimated thetas, which are not individually interpretable. For correlations asymptotic standard errors were obtained by the delta method implemented in TMB. 34 Table 5 Age Case Variation Mean SD 4 1 Per + Trans 222.34 39.45 4 2 Per 222.4 34.22 4 3 Trans 219.5 19.09 15 1 Per + Trans 565.5 88.93 15 2 Per 565.7 87.49 15 3 Trans 570.6 16.59 40 1 Per + Trans 723.2 154.83 40 2 Per 772 115.37 40 3 Trans 714.9 3.46 Table 5: Summary statistics of the frequency distribution of length-at-age for each case and each age (4, 15, and 40). 35 Table 6 Estimate Estimate SE Parameter CI LB CI UB (log) (numerical) (log) Observation Error -7.18 0.0007617 0.21 0.0005005 0.0011592 L∞ deviations -1.39 0.2490753 0.07 0.2165357 0.2865048 (individual) K deviations -1.02 0.3605949 0.08 0.3072787 0.4231621 (individual) L1 deviations -1.51 0.22091 0.05 0.1998876 0.2441433 (individual) L∞ deviations -1.94 0.1437039 0.36 0.0699482 0.2952302 (population) K deviations -1.69 0.1845195 0.4 0.08291 0.4106558 (population) L1 deviations -2.85 0.0578443 0.46 0.0230521 0.1451482 (population) L∞ mean 6.56 706.27169 0.07 614.00311 812.40583 (metapopulation) K mean -2.29 0.1012665 0.09 0.0845849 0.121238 (metapopulation) L1 mean 4.65 104.58499 0.03 98.49443 111.05216 (metapopulation) Process Error -1.45 0.2345703 0.02 0.2253727 0.2441433 L∞ / K correlation -0.81 0.16 (individual) L∞ / L1 correlation -0.03 0.1 (individual) K / L1 correlation 0.28 0.09 (individual) L∞ / K correlation 0.3 0.65 (population) Table 6: Parameter estimates for the male-specific VBGM model. 36 Figure 1: Location of the six populations of lake trout sampled from Lake Superior. 37 Figure 2: Structure of the hierarchical von Bertalanffy growth model. The second level represents Lake Superior-level log-scale means for each parameter. The second level represents population-level deviations, as random effects, for each parameter population (k=6). The third level represents individual-level deviations, as random effects, for each parameter and each individual (n=410). Finally, the individual- and population-level deviations are added to the metapopulation mean to create each individual’s growth parameters; these are then used to predict each individual’s growth trajectory. 38 Figure 3: Process error (transient variation) and observation error in the hierarchical model. Process error is age- and individual-specific, assumed to be normally distributed, and is multiplied by the predicted growth increment. Predicted length-at-age is considered to be true length-at-age. This true length-at-age is subjected to observation error, assumed to be normally distributed, to produce our observed length-at-age from the model. Figure 4: Biphasic growth. Individual growth trajectories that visually displayed what appeared to be a biphasic growth pattern (i.e. linear early growth and asymptotic growth after an inflection point). 39 Figure 5: Median estimates and 95% probability intervals for values of VBGM parameters. Among-population results show the probability intervals of population-specific estimates for L∞, K, and L1 around the meta-population medians. Among-individual results show the probability intervals of individual-specific estimates for L∞, K, and L1 around a selected example population mean (Grand Marais). Probability intervals were calculated using the back-transformed parameter estimate with the error bars equaling 2*standard deviation from that estimate. 40 Figure 6: Frequency distributions of length-at-age for ages 4, 15, and 40. Each simulation case is represented as a curve on the graph: Case 1 = base case (all variation), Case 2 = persistent variation and observation error, Case 3 = transient variation and observation error. 41 Figure 7: Frequency distributions of length-at-age for ages 4, 15, and 40 using the biphasic model. Each simulation case is represented as a curve on the graph: Case 1 = base case (all variation), Case 2 = persistent variation and observation error, Case 3 = transient variation and observation error. 42 Figure 8: Frequency distributions of length-at-age for ages 4, 15, and 40. This is the sex-specific model using only male lake trout. Each simulation case is represented as a curve on the graph: Case 1 = base case (all variation), Case 2 = persistent variation and observation error, Case 3 = transient variation and observation error. 43 LITERATURE CITED 44 LITERATURE CITED Biro, P.A., Post, J.R., 2008. Rapid depletion of genotypes with fast growth and bold personality traits from harvested fish populations. Proc Natl Acad Sci U S A 105, 2919–2922. https://doi.org/10.1073/pnas.0708159105Campana, S.E., 1990. How Reliable are Growth Back-Calculations Based on Otoliths? Canadian Journal of Fisheries and Aquatic Sciences 47, 2219–2227. Caroffino, D.C., Lenart, S.J., 2000. Statistical catch-at-age models used to describe the status of lean lake trout populations in the 1836-Treaty ceded waters of lakes Michigan , Huron , and Superior at the inception of the 2000 Consent Decree A Report Completed by the Modeling Subcommittee. Chavarie, L., Howland, K.L., Harris, L.N., Hansen, M.J., Gallagher, C.P., Harford, W.J., Tonn, W.M., Muir, A.M., Krueger, C.C., 2019. Habitat overlap of juvenile and adult lake trout of Great Bear Lake: Evidence for lack of a predation gradient? Ecology of Freshwater Fish 28, 485–498. https://doi.org/10.1111/eff.12470 Chavarie, L., Muir, A.M., Zimmerman, M.S., Baillie, S.M., Hansen, M.J., Nate, N.A., Yule, D.L., Middel, T., Bentzen, P., Krueger, C.C., 2017. Challenge to the model of lake charr evolution: Shallow and deep-water morphs exist within a small postglacial lake. Biological Journal of the Linnean Society 120, 578–603. https://doi.org/10.1111/bij.12913 DeAngelis, D.L., Rose, K.A., Crowder, L.B., Marschall, E.A., Lika, D., 1993. Fish Cohort Dynamics: Application of Complementary Modeling Approaches. The American Naturalist 142, 604–622. Eveson, J.P., Polacheck, T., Laslett, G.M., 2007. Consequences of assuming an incorrect error structure in von Bertalanffy growth models: a simulation study. Canadian Journal of Fisheries and Aquatic Sciences 64, 602–617. Flinn, S.A., Midway, S.R., 2021. Trends in growth modeling in fisheries science. Fishes 6, 1–18. https://doi.org/10.3390/fishes6010001 Hansen, M.J., 1999. Lake trout in the Great Lakes: Basin-wide stock collapse and binational restoration, in: Taylor, W.W., Paola Ferreri, C. (Eds.), Great Lakes Fishery Policy and Management: A Binational Perspective. Michigan State University Pres, East Lansing, MI, pp. 417–453. Hansen, M.J., Guy, C.S., Bronte, C.R., Nate, N.A., 2021. Life History and Population Dynamics, in: Muir, A.M., Krueger, C.C., Hansen, M.J., Riley, S.C. (Eds.), The Lake Charr Salvelinus Namaycush: Biology, Ecology, Distribution, and Management. Fish & Fisheries. 45 Hansen, M.J., Madenjian, C.P., Selgeby, J.H., Helser, T.E., 1997. Gillnet selectivity for lake trout (Salvelinus namaycush ) in Lake Superior. Canadian Journal of Fisheries and Aquatic Sciences 54, 2483–2490. https://doi.org/10.1139/f97-156 Hansen, M.J., Nate, N.A., Chavarie, L., Muir, A.M., Zimmerman, M.S., Krueger, C.C., 2016a. Life history differences between fat and lean morphs of lake charr (Salvelinus namaycush) in Great Slave Lake, Northwest Territories, Canada. Hydrobiologia 783, 21–35. https://doi.org/10.1007/s10750-015-2633-2 Hansen, M.J., Nate, N.A., Krueger, C.C., Zimmerman, M.S., Kruckman, H.G., Taylor, W.W., 2012. Age, growth, survival, and maturity of lake trout morphotypes in lake mistassini, quebec. Trans Am Fish Soc 141, 1492–1503. https://doi.org/10.1080/00028487.2012.711263 Hansen, M.J., Nate, N.A., Muir, A.M., Bronte, C.R., Zimmerman, M.S., Krueger, C.C., 2016b. Life history variation among four lake trout morphs at Isle Royale, Lake Superior. Journal of Great Lakes Research 42, 421–432. https://doi.org/10.1016/j.jglr.2015.12.011 Hart, D.R., 2001. Individual-based yield-per-recruit analysis, with an application to the Atlantic sea scallop, Placopecten magellanicus. Canadian Journal of Fisheries and Aquatic Sciences 58, 2351–2358. https://doi.org/10.1139/cjfas-58-12-2351 Hart, D.R., Chute, A.S., 2009. Estimating von Bertalanffy growth parameters from growth increment data using a linear mixed-effects model, with an application to the sea scallop Placopecten magellanicus. ICES Journal of Marine Science 66, 2165–2175. https://doi.org/10.1093/icesjms/fsp188 He, J.X., Bence, J.R., Madenjian, C.P., Pothoven, S.A., Dobiesz, N.E., Fielder, D.G., Johnson, J.E., Ebener, M.P., Adam Cottrill, R., Mohr, L.C., Koproski, S.R., 2015. Coupling age- structured stock assessment and fish bioenergetics models: A system of time-varying models for quantifying piscivory patterns during the rapid trophic shift in the main basin of lake huron. Canadian Journal of Fisheries and Aquatic Sciences 72, 7–23. https://doi.org/10.1139/cjfas-2014-0161 Hixon, M.A., Johnson, D.W., Sogard, S.M., 2014. BOFFFFs: On the importance of conserving old-growth age structure in fishery populations. ICES Journal of Marine Science 71, 2171– 2185. https://doi.org/10.1093/icesjms/fst200 James, I.R., 1991. Estimation of von Bertalanffy Growth Curve Parameters from Recapture Data. Biometrics 47, 1519–1530. Kimura, D.K., Shimada, A.M., Lowe, S.A., 1993. Estimating von Bertalanffy growth parameters of sablefish Anoplopoma fimbria and Pacific cod Gadus macrocephalus using tag-recapture data. Fishery Bulletin 91, 271–280. 46 Kirkwood, G.P., Somers, I.F., 1984. Growth of two species of tiger prawn, penaeus esculentus and p. Semisulcatus, in the western gulf of carpentaria. Marine and Freshwater Research 35, 703–712. https://doi.org/10.1071/MF9840703 Kraak, S.B.M., Haase, S., Minto, C., Santos, J., Anderson, E., 2019. The Rosa Lee phenomenon and its consequences for fisheries advice on changes in fishing mortality or gear selectivity. ICES Journal of Marine Science 76, 2179–2192. https://doi.org/10.1093/icesjms/fsz107 Kristensen, K., Bell, B., Skaug, H., Magnusson, A., Berg, C., Nielsen, A., Maechler, M., Michelot, T., Brooks, M., Forrence, A., Moesgaard, C., Monnahan, C., 2022. Template Model Builder: A General Random Effect Tool Inspired by “ADMB.” Kristensen, K., Nielsen, A., Berg, C.W., Skaug, H., Bell, B.M., 2016. TMB: Automatic differentiation and laplace approximation. Journal of Statistical Software 70. https://doi.org/10.18637/jss.v070.i05 Kristiansen, T.S., Svåsand, T., 1998. Effect of size-selective mortality on growth of coastal cod illustrated by tagging data and an individual-based growth and mortality model. Journal of Fish Biology 52, 688–705. https://doi.org/10.1006/jfbi.1997.0611 Lenart, S.J., Caroffino, D.C., 2019. Technical Fisheries Committee Administrative Report 2019: Status of Lake Trout and Lake Whitefish Populations in the 1836 Treaty-Ceded Waters of Lakes Superior, Huron and Michigan, with Recommended Yield and Effort Levels for 2019. Lester, N.P., Shuter, B.J., Abrams, P.A., 2004. Interpreting the von Bertalanffy model of somatic growth in fishes: The cost of reproduction. Proceedings of the Royal Society B: Biological Sciences 271, 1625–1631. https://doi.org/10.1098/rspb.2004.2778 Lester, N.P., Shuter, B.J., Venturelli, P., Nadeau, D., 2014. Life-history plasticity and sustainable exploitation: A theory of growth compensation applied to walleye management. Ecological Applications 24, 38–54. https://doi.org/10.1890/12-2020.1 Martínez-Garmendia, J., 1998. Simulation analysis of evolutionary response of fish populations to size-selective harvesting with the use of an individual-based model. Ecological Modelling 111, 37–60. https://doi.org/10.1016/S0304-3800(98)00093-3 Morrongiello, J.R., Thresher, R.E., Hobbs, N.T., 2015. A statistical framework to explore ontogenetic Growth variation among individuals and populations: A marine fish example. Ecological Monographs 85, 93–115. https://doi.org/10.1890/13-2355.1 Muir, A.M., Bronte, C.R., Zimmerman, M.S., Quinlan, H.R., Glase, J.D., Krueger, C.C., 2014. Ecomorphological Diversity of Lake Trout at Isle Royale, Lake Superior. Trans Am Fish Soc 143, 972–987. https://doi.org/10.1080/00028487.2014.900823 Muir, A.M., Krueger, C.C., Hansen, M.J., 2012. Re-Establishing Lake Trout in the Laurentian Great Lakes: Past, Present, and Future, in: Taylor, W.W. (Ed.), Great Lakes Fisheries Policy 47 and Management: A Binational Perspective (2ed). Michigan State University Press, pp. 533– 588. Muir, A.M., Krueger, C.C., Hansen, M.J., Riley, S.C. (Eds.), 2021. The Lake Charr Salvelinus namaycush: Biology, Ecology, Distribution, and Management, Fish & Fisheries Series. Fish & Fisheries. Ogle, D.H., Brenden, T.O., McCormick, J., 2012. Growth Estimation : Growth Models and Statistical Inference, in: Quist, M., Isermann, D. (Eds.), Age and Growth of Fishes: Principles and Techniques. American Fisheries Society, Bethesda, Maryland, pp. 265–359. Ortiz de Zarate, V., Babcock, E.A., 2015. Estimating individual growth variability in albacore (Thunnus alalunga) from the North Atlantic stock; aging for assessment purposes. Parma, A.M., Deriso, R.B., 1990. Dynamics of Age and Size Composition in a Population Subject to Size-Selective Mortality: Effects of Phenotypic Variability in Growth. Canada Journal of Fisheries and Aquatic Sciences 47. Pauly, D., 1980. On the interrelationships between natural mortality, growth parameters, and mean environmental temperature in 175 fish stocks. J. Cons. int. Explor. Mer 39, 175–192. Peacor, S.D., Bence, J.R., Pfister, C.A., 2007. The effect of size-dependent growth and environmental factors on animal size variability. Theoretical Population Biology 71, 80–94. https://doi.org/10.1016/j.tpb.2006.08.005 Pfister, C.A., Stevens, F.R., 2002. The Genesis of Size Variability in Plants and Animals. Ecology 83, 59. https://doi.org/10.2307/2680121 Pilling, G.M., Kirkwood, G.P., Walker, S.G., 2002. An improved method for estimating individual growth variability in fish, and the correlation between von Bertalanffy growth parameters. Canadian Journal of Fisheries and Aquatic Sciences 59, 424–432. https://doi.org/10.1139/f02-022 Punt, A.E., Allen Akselrud, C., Cronin-Fine, L., 2017. The effects of applying mis-specified age- and size-structured models. Fisheries Research 188, 58–73. https://doi.org/10.1016/j.fishres.2016.11.017 Punt, A.E., Hobday, D., Gerhard, J., Troynikov, V.S., 2006. Modelling growth of rock lobsters, Jasus edwardsii, off Victoria, Australia using models that allow for individual variation in growth parameters. Fisheries Research 82, 119–130. https://doi.org/10.1016/j.fishres.2006.08.003 Quince, C., Abrams, P.A., Shuter, B.J., Lester, N.P., 2008a. Biphasic growth in fish I: Theoretical foundations. Journal of Theoretical Biology 254, 197–206. https://doi.org/10.1016/j.jtbi.2008.05.029 48 Quince, C., Abrams, P.A., Shuter, B.J., Lester, N.P., 2008b. Biphasic growth in fish II: Empirical assessment. Journal of Theoretical Biology 254, 197–206. https://doi.org/10.1016/j.jtbi.2008.05.029 Quinn, T.J., Deriso, R.B., 1999. Quantitative Fish Dynamics. Oxford University Press, New York. Ricker, W.E., 1969. Effects of Size-Selective Mortality and Sampling Bias on Estimates of Growth, Mortality, Production, and Yield. Journal of the Fisheries Research Board of Canada 26. Sainsbury, K.J., 1980. Effect of Individual Variability on the von Bertalanffy Growth Equation. Canadian Journal of Fisheries and Aquatic Sciences 37, 241–247. Stawitz, C.C., Essington, T.E., Branch, T.A., Haltuch, M.A., Hollowed, A.B., Spencer, P.D., 2015. A state-space approach for detecting growth variation and application to North Pacific groundfish. Canadian Journal of Fisheries and Aquatic Sciences 72, 1316–1328. https://doi.org/10.1139/cjfas-2014-0558 Stawitz, C.C., Haltuch, M.A., Johnson, K.F., 2019. How does growth misspecification affect management advice derived from an integrated fisheries stock assessment model? Fisheries Research 213, 12–21. https://doi.org/10.1016/j.fishres.2019.01.004 Thorson, J.T., Minte-Vera, C. v., 2016. Relative magnitude of cohort, age, and year effects on size at age of exploited marine fishes. Fisheries Research 180, 45–53. https://doi.org/10.1016/j.fishres.2014.11.016 Tsehaye, I., Jones, M.L., Bence, J.R., Brenden, T.O., Madenjian, C.P., Warner, D.M., 2014. A multispecies statistical age-structured model to assess predator-prey balance: Application to an intensively managed Lake Michigan pelagic fish community. Canadian Journal of Fisheries and Aquatic Sciences 71, 627–644. https://doi.org/10.1139/cjfas-2013-0313 Vigliola, L., Meekan, M.G., 2009. The Back-Calculation of Fish Growth From Otoliths, Media. https://doi.org/10.1007/978-1-4020-5775-5 Vincenzi, S., Mangel, M., Crivelli, A.J., Munch, S., Skaug, H.J., 2014. Determining Individual Variation in Growth and Its Implication for Life-History and Population Processes Using the Empirical Bayes Method. PLoS Computational Biology 10. https://doi.org/10.1371/journal.pcbi.1003828 Wang, Y.-G., 1998. An improved Fabens method for estimation of growth parameters in the von Bertalanffy model with individual asymptotes. Canadian Journal of Fisheries and Aquatic Sciences 55, 397–400. https://doi.org/10.1139/cjfas-55-2-397 49 Webber, D.N., Thorson, J.T., 2016. Variation in growth among individuals and over time: A case study and simulation experiment involving tagged Antarctic toothfish. Fisheries Research 180, 67–76. https://doi.org/10.1016/j.fishres.2015.08.016 Weisberg, S., Spangler, G., Richmond, L.S., 2010. Mixed effects models for fish growth. Canadian Journal of Fisheries and Aquatic Sciences 67, 269–277. https://doi.org/10.1139/F09-181 Wellband, K., Baillie, S.M., Bentzen, P., Bernatchez, L., 2021. Genetic Diversity, in: Muir, A.M., Krueger, C.C., Hansen, M.J., Riley, S.C. (Eds.), The Lake Charr Salvelinus Namaycush: Biology, Ecology, Distribution, and Management. Fish & Fisheries. Wolf, M., Weissing, F.J., 2012. Animal personalities: Consequences for ecology and evolution. Trends in Ecology and Evolution. https://doi.org/10.1016/j.tree.2012.05.001 Xiao, Y., 1994. von Bertalanffy growth models with variability in, and correlation between, K and Linfinity. Canadian Journal of Fisheries and Aquatic Sciences 51, 1585–1590. https://doi.org/10.1139/f94-157 Zhang, Z., Lessard, J., Campbell, A., 2009. Use of Bayesian hierarchical models to estimate northern abalone, Haliotis kamtschatkana, growth parameters from tag-recapture data. Fisheries Research 95, 289–295. https://doi.org/10.1016/j.fishres.2008.09.035 50 CHAPTER TWO: EFFECTS OF INDIVIDUAL VARIABILITY IN GROWTH ON MANAGEMENT REFERENCE POINTS: AN APPLICATION FOR LAKE TROUT (SALVELINUS NAMAYCUSH) IN LAKE SUPERIOR AUTHORS Elizabeth Stebbins1 James R. Bence1 Travis O. Brenden1 1 Quantitative Fisheries Center, Michigan State University ABSTRACT Individual variability in fish growth is a well-documented phenomenon and several methods have been developed to understand its influence on estimates of population parameters in fisheries, showing that aggregating growth data by age can cause bias or inaccuracies in estimates. We developed a simulation framework to assess the effect of individual variability in growth on basic per-recruit models used in management. We applied our framework based on estimated parameters for a typical population of lake trout (Salvelinus namaycush) in Lake Superior. We use two approaches to calculate reference points: a ‘standard’ approach that assumes access to good data but ignores individual variability, and a ‘true’ approach that accounts for the interaction between size-selective mortality and growth variability for individuals. We show that the standard method consistently overestimates yield-per-recruit, especially at higher levels of instantaneous fishing mortality, compared to the true method. Further, the standard method underestimates spawning stock biomass-per-recruit at low levels of fishing intensity. Because of the differential biases as a function of fishing mortality, reference points based on the standard approach showed anti-conservative bias (i.e., suggesting a higher level of F could be sustained). Therefore, ignoring growth variation could result in 51 mismanagement of stocks or risk. This study contributes to a body of work that attempts to understand the implications of individual variability in growth to management reference points and proposes a potential solution framework in the form of an individual-specific model. 52 INTRODUCTION Quantifying fish growth though the fitting of fish growth models to size-at-age or length increment data is an important aspect of fishery science. Estimating growth parameters can inform population dynamics like mortality rates and consequently can be informative when it comes to setting management policies. Common growth models like the von Bertalanffy growth model (VBGM) and its various parameterizations can predict incremental growth (in length or in weight) over time and have been applied to a wide variety of fish species. In fishery science, the VBGM is by far the most common growth model and often is used to estimate population-level parameters, although it can also be used to model individual-specific growth so as long as informative data are available (Flinn & Midway, 2021; Sainsbury, 1980). Estimating growth processes correctly is important because many of the biological reference points used in analyses like statistical catch-at-age models or per-recruit models require information about the weight at age of individuals throughout their lives (Quinn & Deriso, 1999a). Analyses that calculate the yield or reproductive potential of fish populations are commonly used to assess harvest potential of a population. These analyses determine the amount of harvest produced by an individual fish (yield-per-recruit) or the number of eggs or spawning stock produced by an individual fish (spawning stock biomass per-recruit) and can form the basis for calculating biological reference points that inform harvest regulations. A basic yield-per-recruit (YPR) or spawning stock biomass-per-recruit (SSBR) analysis incorporates information about growth (usually as mean weight-at-age of a harvested fish or mean weight-at-age at time of spawning), instantaneous fishing mortality F, size-selectivity of the gear used, natural mortality M, and probability of a fish being mature. With access to these data, analysts can calculate how yield or spawning biomass per recruitment will change with fishing, and metrics that implicitly 53 make assumptions about how recruitment will change, that inform on how the level of F may affect population demography, such as abundance and size structure. Basic YPR and SSBR analyses generally aggregate individual values into population means (for example, mean weight-at-age) (Quinn & Deriso, 1999a). Although aggregating individuals simplifies analyses, this also ignores the substantial variability that is frequently found in fish populations like that estimated in Chapter 1. According to Parma and Deriso (1990), per-recruit estimates can be biased if analyses do not account for the effects of individual variability and size-selective mortality on growth, though their analysis assumed the fishing season was short enough such that growth and natural mortality did not occur alongside fishing mortality, and based their model on parameter estimates from Pacific halibut (Hippoglossus stenolepis) (Parma & Deriso, 1990). Further, size-selective fishing or other sources of size-selective mortality can alter the size-structure of a stock by preferentially removing fast-growing, larger fish, resulting in slow-growing cohorts dominating the stock at higher ages. This result is termed the ‘Rosa Lee phenomenon’ (Lee, 1912) whereby back-calculated lengths at ages tend to be less than the length at age from fish captured at that age. The effects of this on estimation of population parameters is generally thought to be significant (Kraak et al., 2019b; Ricker, 1969b). While the Rosa Lee phenomenon is not the focus of our study, it encapsulates the interaction between size-selective fishing mortality and individual variability in growth and how it might influence management policies for a population. A key input to basic YPR and SSBR analyses is the mean weight-at-age of fish in a stock. For our purposes, we focus on the discrete-time generic, or Bell-Thompson, method outlined in Quinn & Deriso (1999), though we acknowledge that there are more complex per-recruit models that can be modified to incorporate the effects of increased fishing intensity or other parameters 54 and we discuss this later (Quinn & Deriso, 1999b). By aggregating weight or length of individuals into average values for each age, an implicit assumption is made the average size at age represents growth of all individuals (i.e., that every individual has the same growth trajectory). The potential problem with this approach is that fish often show different growth trajectories within a population. Logically, a size-selective mortality factor (for example, gill nets that exclude fish below a certain size) would remove some individuals sooner than others depending on how fast the individual reaches the size of entry to the fishery, thus changing the distribution of growth characteristics of the survivors. When length-at-age is aggregated in a model, this nuance, and any subsequent effects on the probability of an individual fish being harvested or surviving, is not explicitly modeled. There have been several models developed to assess the interactions between individual variability, size-selective mortality, and per-recruit model calculations. Many of these models assume that this variability is due to intrinsic characteristics of the fish (persistent growth) as opposed to environmental conditions (transient growth) or a combination of the two (Kristiansen & Svåsand, 1998; Lowerre-Barbieri et al., 1998; Webber & Thorson, 2016). Parma & Deriso (1990) discuss the difficulty of disentangling the relative contribution of the two sources of variability, and our previous study presents one method for attempting to estimate them (Chapter 1). Persistent growth can potentially interact significantly with fishing to alter size structure of populations, genetic composition, or evolutionary trajectories of a stock (Kuparinen et al., 2009; Wolf & Weissing, 2012). Relatedly, the effects of size-selective mortality and persistent growth variability on evolution of a fish stock have been studied but are generally not well understood (Martínez-Garmendia, 1998). In some cases, yield has been shown to be underestimated when individual variability is ignored (Kristiansen & Svåsand, 1998; Parma & Deriso, 1990), but is highly 55 dependent on fishing intensity. Estimates of spawning stock at a level of exploitation relative to unfished spawning stock have been shown to be overestimated when individual variability is ignored (Stawitz et al., 2019). Therefore, it is also important to consider how exploited the stock under question is. Clearly, incorporating individual variability in growth when calculating biological reference points and yield- or spawning stock biomass models is important to consider explicitly, and is becoming easier with advances in modeling and increases in computational power. In a previous study, we developed a mixed-effects, hierarchical model to estimate the magnitude and characteristics of individual- and population-level variability in von Bertalanffy growth model parameters L∞ (asymptotic length), K (Brody growth coefficient, or how fast the fish reaches its asymptotic length), and L1 (length at age 1). Our analysis showed that individual- and population-level variability existed for these parameters and for some parameters this variability could be substantial (Chapter 1). Here we present an analysis that uses the parameter estimates obtained in the previous study to simulate the growth of individuals and, paired with size- selective mortality and maturity patterns, their population dynamics. We then calculate observed values of YPR, SSBR, and biological reference points using a ‘standard’ approach based on mean size at age and other inputs obtained under a reference level of fishing and compare this with true values that take into account variation in growth among individuals and how this interacts with size selective fishing. We outline our approach, which uses different growth morphs to represent individual variability in growth (in a method similar to (Martínez- Garmendia, 1998), below. We hypothesized that the standard approach using averaged sizes at age would differ from the true values, and that the magnitude of F would be an important determinant of the magnitude of difference in these values. Ultimately, our goal was to expand 56 knowledge about how size-selective mortality interacts with individual variability in growth and how this may be of interest to managers when setting harvest policies. 57 METHODS Overview We assessed the influence of growth variability, here represented by ‘growth morphs’ and their different growth trajectories, on per-recruit analyses. We allowed for variation in growth to interact with length-based selectivity and maturity functions, the former of which preferentially removes some growth morphs from the population. The basic idea underlying growth morphs is each morph represents the persistent characteristics associated with growth of a larger number of individuals with similar VBGM parameters, and our analysis accounted for how common the parameters associated with each growth morph was expected to be. We simulated the growth of one individual (super-individual) for each growth morph. While different individuals from the same growth morph will vary in their growth due to transient variation in growth, we showed in Chapter 1 that such transient growth was relatively unimportant in producing variation in size at age for lake trout in Lake Superior and thus believe using one super-individual to represent a growth morph is reasonable. In support of our decision to simulate only one individual we found very similar results with transient variation set to zero, showing that transient variation was having little influence However, if transient variation were expected to be more dominant, we would have considered modeling several fish per node. We simulated growth trajectories for 1000 fish, each representing a growth morph, using ten values of each of the VBGM parameters L∞, K and L1 in unique combinations (described below). We then grew each individual from age 1 to age 50 in daily growth increments and created individual-specific length-based selectivity and maturity schedules modeled using parameterizations for Lake Superior lake trout (Figure 9). This produced individual-specific time series in daily increments for length, weight, maturity, and selectivity. We combined these time 58 series with fully selected fishing mortality F and natural mortality M to identify the relative frequencies of each super-individual that survived to each day and used that to calculate the probability of each super-individual’s survival and probability of being harvested on a given day. With these data representing the “truth”, we then performed per-recruit analyses using a ‘standard’ approach requiring average weight-at-age, proportion mature-at-age, and selectivity- at-age, each in the form of annual values. To obtain these values, we assumed the analyst had good knowledge about the system and calculated weighted averages using each growth morph’s relative frequency expected by the multivariate normal distribution estimated in Chapter 1. We then analyzed the performance of the standard approach by comparing the per-recruit values and reference points calculated from them versus the true values that accounted for among individual variability. Parameters, their symbols, and their descriptions are explained below and can be referenced in Table 7. For all analyses we used R statistical software, version 4.0.2 (R Core Team, 2020). Individual variation in growth We assumed that individual fish grew according to a VBGM where the expected growth increment any given day was determined by the parameters, transient growth variation, and the starting length of a fish that day. We used the standard VBGM parameters L∞ and K to represent asymptotic size and the Brody growth coefficient, but instead of t0, we used L1 (length of an age-1 fish at the beginning of day 1) because it is more biologically interpretable (Eq. 1). 𝐿(1, 𝑖 ) = 𝐿1 𝑖 (1) Persistent variability was simulated by allowing each growth morph to have a unique combination of growth parameters. Our approach was to construct a three-dimensional grid of possible growth parameters (10 L∞ values x 10 K values x 10 L1 values), so 1000 growth 59 morphs. We assigned one individual fish (super-indvidual) to each growth morph. In the previous chapter we fit a hierarchical growth model describing the trivariate distribution of the VBGM parameters among individuals and obtained estimates of transient variation represented by a process error (Chapter 1). The parameter estimates describing the trivariate distribution was used at a later stage to weight the fish representing each growth morph to account for the expected differences in the frequencies of the different growth morphs. Depending on the estimated distribution of each VBGM parameter from the previous chapter (log-scale mean and standard deviation), we calculated the 10 values for each parameter by calculating the different quantiles given the other parameters were at their mean values. We used the hyper-population mean over populations as being representative of a typical population and the standard deviation for among-individual variation within a population as defining a typical distribution for parameters. This was to ensure we evaluated a wide range of potential values encompassing the plausible growth morphs. The specific approach, which has some complexity, was merely one method to define a repeatable method for choosing values that spanned the range of plausible values. In particular, we used the mean and standard deviation for each parameter (on a log scale) to calculate 10 bins of equal probability for a parameter (conditional on the other log scale growth parameters being at their average values). More specifically, we calculated the quantile values for probabilities 0.1 – 0.9, and these defined the 10 bins (with the lower bound of the first bin and upper bound of the last bin based on the quantiles for 0.01% and 0.999%). We then set for each growth parameter and bin an associated parameter value that was the mean for that parameter for that bin, based on the normal distribution for that log scale parameter (with other parameters at their mean values). 60 The log scale mean for each bin for each log scale parameter was calculated based on known results for truncated normal distributions (Burkardt, 2014) by: 𝜙 (𝛼 ) − 𝜙 (𝛽 ) (2.1) 𝑚𝑝 = 𝜇 + 𝜎 𝛷 (𝛽) − 𝛷(𝛼) with 𝛼 and 𝛽 calculated based on the lower (a) and upper (b) limits of the interval (i.e., the bounds calculated previously) (Eq. 2.1). 𝑎−𝜇 𝛼= (2.2) 𝜎 𝑏 −𝜇 (2.3) 𝛽= 𝜎 and μ and σ the estimated mean and standard deviation for the log-scale parameters describing among individual variation in each parameter (taken from the estimation model). Here 𝜙(𝛼 ) and 𝜙(𝛽) are the probability densities of 𝛼 and 𝛽 given a standard normal distribution (mean of 0 and standard deviation of 1), and 𝛷(𝛽) and 𝛷(𝛼) are the cumulative densities of 𝛼 and 𝛽 given a standard normal distribution. Simulated growth We simulated the growth of an individual fish I using one of the 1000 unique combinations of log scale VBGM parameters defined above (growth morphs). These morph-specific growth parameters influence the trajectory of a fish’s size-at-age over the course of its life, representing intrinsic persistent variation in growth. We started growth at age-1 and then calculated the expected growth increment for a day d using the Fabens parameterization of the VBGM (Fabens, 1965) (Eq. 3.1). Individual fish from the same growth morph (i.e., with the same VBGM parameters) will follow somewhat different growth trajectories because of transient process growth variation. As explained previously we simulated just one super-individual per growth morph. To model transient growth, which is external and year-to-year, we simulated process error as an individual- and year-specific error multiplied by the expected growth increment (Eq. 61 3.2). Because the error term is year-specific but applied to a daily increment, we divided the error term by 365 which assumes that the transient variability in growth occurs equally over all days of the year. 𝛿 (𝑎,𝑖) (3.1) ΔL(𝑑, 𝑖 ) = (𝐿∞𝑖 − 𝐿(𝑑 − 1, 𝑖)(1 − exp(−𝐾𝑖 )) ∗ 365 The process error is modeled on the log scale as normally distributed with a bias correction: −𝜎 2 (3.2) log(𝛿 (𝑎, 𝑖 )) ~𝑁 ( 𝛿 , 𝜎𝛿2 ) 2 For day 1 of an age-1 fish, we set an individual’s length to 𝐿1 𝑖 (Eq. 1). For all days after day 1 (d > 1), length-at-day was calculated by adding the increment to the length at the previous day: 𝐿 (𝑑, 𝑖 ) = 𝐿 (𝑑 − 1, 𝑖 ) + ΔL(𝑑, 𝑖 ), 𝑑 > 1 (3.3) To calculate weight as a function of length, we used an allometric relationship: 𝑤 = 𝜆𝐿𝛽 (3.4) We obtained estimates of the weight-length parameters by regressing natural log weight (g) on log total length (mm) at time of capture for the same fish that provided the data used to fit the hierarchical model. This produced estimates of β (slope of regression) and 𝜆 (log intercept of regression) equal to 3.19 and 2.35e-06, respectively. We then compared these values to those calculated for Lake Huron (He et al., 2008) and Lake Michigan (Matthias, Bence, & Clark in prep) trout to ensure they were reasonable. For simplicity, we fixed β and 𝜆 for our simulations, but recognize that this has a significant impact on later yield and spawning stock biomass calculations. 62 Calculating weights (𝛾) for each growth morph After growing the 1000 growth morphs to age 50 we needed to calculate the probability of each growth morph’s existence in a realistic population. Because the 10 different values of each VBGM parameter were not equally likely to occur in a population (i.e., combination of some of the more extreme values were less likely to occur than those closer to the mean), we calculated the probability of observing each of the 1000 combinations of parameters. These were then used in later analyses to take the weighted average of metrics. Simulated selectivity and maturity schedules We modeled selectivity patterns and maturity schedules as simple logistic functions of length. For the lake trout gillnet fishery, a double-logistic selectivity function is often used (Lenart & Caroffino, 2019; Sitar & He, 2006). We used a simple logistic curve to induce simplicity in the model, and note that the selectivity function in this framework could easily be adapted to the specific pattern and gear type of the fishery (Kuparinen et al., 2009). For selectivity: 1 (4) 𝑠𝑒𝑙𝐿,𝑑 = (1 + exp(−𝛼(𝐿𝑑 − 𝐿𝐸 ))) 𝐿𝐸 is the length at age-of-entry to the fishery, which we set at 528 mm. For 𝛼, we used a value of 0.02 (Bence et al., 2003). We applied the selectivity function to the length at each day over a 50- year life span. The logistic maturity schedule was modeled as: 1 (5) 𝑝𝑚𝑎𝑡𝐿,𝑑 = (1 + exp(−𝛽(𝐿𝑑 − 𝐿50 ))) For 𝐿50, we used a value of 558.08 mm and chose a value of 0.02 for 𝛽 which were obtained from the 1836 treaty waters Modeling Subcommittee stock assessment models for Lake Superior 63 length-based maturity schedule to reflect actual population models used for the system (Shawn Sitar, Michigan DNR, personal communication). This function was applied to the same length at each day for a fish that lived to 50 years as was done for selectivity. Calculations of age-specific averages from known data To assess the accuracy of a standard approach to individual variability represented by growth morphs, we made assumptions about the quality and type of data to which an analyst had access. For example, to calculate yield-per-recruit and spawning stock biomass-per-recruit using methods outlined in the discrete-time generic method, the analyst would need selectivity-at-age, proportion mature-at-age, average weight-at-age at time of harvest, and average weight-at-age at time of spawning (Quinn & Deriso, 1999b). We used a series of calculations to produce these age- based patterns from known simulated values for selectivity and maturity pattern in daily increments and daily length and weight increments. Total mortality, Z, was calculated using the selectivity-at-day pattern, a baseline annual fishing mortality F of 0.1 (that applies when selectivity is 1), and a baseline annual natural mortality M 1 of 0.1 (that applied across all days). We multiplied total annual mortality Z by to convert to a 365 daily rate: 1 (6.1) 𝑍𝑑,𝑖 = (𝑠𝑒𝑙𝑑,𝑖 ∗ 𝐹 + 𝑀) 365 First, abundance per recruit at the beginning of each day was computed for each super- individual. Recruits (abundance at the start of day 1) was set equal to 1, which means that abundance in each subsequent day was the abundance per initial recruit at the start of day 1. 𝑁𝑑+1,𝑖 = 𝑁𝑑,𝑖 ∗ 𝑒 (−𝑍𝑑,𝑖 ) (6.2) 64 This is equivalent to the probability of a fish being alive at the beginning of a day d. The probability of each super-individual I being harvested on a given day d was then computed as: 𝐹𝑑,𝑖 ∗ (1 − 𝑒 −𝑍𝑑,𝑖 ) (6.3) 𝑃(𝐻)𝑑,𝑖 = ∗ 𝑁𝑑,𝑖 𝑍𝑑,𝑖 These values were used to calculate average values of abundance (𝑁𝑎 ), catch (𝐶𝑎 ), selectivity (𝑠𝑒𝑙𝑎 ), and proportion mature (𝑝𝑚𝑎𝑡𝑎 ) at age a. Abundance at age a, 𝑁𝑎 , in years, was the abundance at the start of the first day of the year for that age. Thus, 𝑁𝑎 was calculated as the weighted average over growth morphs: 𝐼 (6.4) 𝑁𝑎 = ∑ 𝛾𝑖 ∗ 𝑃(𝐴)𝑎,𝑖 𝑖 Where 𝑃 (𝐴)𝑎,𝑖 = 𝑁(𝑎−1)∗365+1 and I = the total number of growth morphs (I = 1000). To calculate the probability of observing a combination of VBGM parameters (𝛾), the estimated trivariate normal distribution for the VBGM parameter was used Chapter 1. Weights were the probability densities of each combination of {L∞, K, and L1} with multivariate normal means, standard deviations, and covariances, normalized to sum to 1.0 over all 1000 growth morphs. This effectively used a discretized version of the estimated distribution of growth parameters to describe the probability a given individual fish, in reality, would be from each of the 1000 growth morphs. Probability densities were calculated used the dmvnorm function from the mvtnorm package in R (Genz et al., 2020). The total catch of age-a fish during a year for fish alive at the start of the first day of the year was computed as: 𝐼 𝑑=365(𝑎−1)+365 (6.5) 𝐶𝑎 = ∑ 𝛾𝑖 ∗ ∑ 𝑃(𝐻)𝑑,𝑖 𝑖 𝑑=365(𝑎−1)+1 65 We assumed the analyst had access to good 𝑁𝑎 and 𝐶𝑎 values, and therefore could infer a selectivity-at-age 𝑠𝑒𝑙𝑎 using the Baranov catch equation: 𝐹 ∗ 𝑠𝑒𝑙𝑎 (6.6) 𝐶𝑎 = (1 − 𝑒 −𝐹∗ 𝑠𝑒𝑙𝑎 −𝑀 )𝑁𝑎 𝐹 ∗ 𝑠𝑒𝑙𝑎 + 𝑀 By moving 𝐶𝑎 to the right side of this equation and solving for the root using an iterative search, the analyst would be able to obtain 𝑠𝑒𝑙𝑎 consistent with values of 𝐶𝑎 and 𝑁𝑎 . The uniroot function from the stats package in R was used (R Core Team, 2020). We are not implying that actual analysts will generally obtain selectivity patterns in this way but rather that the best an analyst could do is have a good understanding of survival and catch produced per recruit and calculate selectivity by year of age to be consistent with that information. The final information needed for the analyst to implement standard YPR and SSBR analyses were average weight-at-age at harvest (𝑊𝑎𝐻 ) and average weight-at-age at spawning (𝑊𝑎𝑠𝑝 ). Weight-at-day was estimated for each growth morph (𝑊𝑑,𝑖 ) by applying Eq. 3.4 to length increments calculated for each growth morph. First, growth during a day was approximated by assuming geometric growth and average weight-at-day for a growth morph (𝑊 ̅𝑑,𝑖 ) as weight at midday: 𝑊̅𝑑,𝑖 = 𝑊𝑑,𝑖 ∗ 𝑒 𝐺𝑑,𝑖 ∗0.5 (7.1) Where 𝐺𝑑,𝑖 is: 𝑊𝑑+1 (7.2) 𝐺𝑑,𝑖 = ln ( ) 𝑊𝑑 Using 𝑊 ̅ 𝑑,𝑖 , we calculated 𝑊𝑎 using the following equation: 𝐻 (𝑎−1)+365 ∑𝑑=365 𝑑=365(𝑎−1)+1 𝑃(𝐻)𝑑,𝑖 ∗ 𝑊𝑑,𝑖 ̅̅̅̅̅ (7.3) ∑𝐼𝑖 𝛾𝑖 ∗ 𝑃(𝐻 )𝑎,𝑖 ∗ (𝑎−1)+365 ∑𝑑=365 𝑑=365(𝑎−1)+1 𝑃(𝐻)𝑑,𝑖 𝑊𝑎𝐻 = ∑𝐼𝑗 𝛾𝑗 ∗ 𝑃(𝐻)𝑎,𝑖 Where: 66 𝑑=365(𝑎−1)+365 (7.4) 𝑃(𝐻)𝑎 = ∑ 𝑃(𝐻)𝑑 𝑑=365(𝑎−1)+1 This first calculates the average harvested weight (in g) of a growth morph for each age, weighted in proportion to the probability of harvest within the year. Average weight-at-age at spawning (𝑊𝑎𝑠𝑝 ) was calculated using values from the day of spawning (day 304 of each year, so sp = 304, or October 31st): ∑𝐼𝑖 𝛾𝑖 ∗ 𝑃(𝐴)365(𝑎−1)+𝑠𝑝,𝑖 ∗ 𝑊365(𝑎−1)+𝑠𝑝,𝑖 (7.5) 𝑊𝑎𝑠𝑝 = 𝐼 ∑𝑗 𝛾𝑗 ∗ 𝑃(𝐴)365(𝑎−1)+𝑠𝑝,𝑖 The day of spawning was the same as the 1836 treaty waters Modeling Subcommittee stock assessment models for Lake Superior (Shawn Sitar, Michigan DNR, personal communication). For probability of being alive on the day of spawning within a year, weight at the start of the day on day of spawning was used, rather than weight at the middle of the day. Finally, average proportion-mature-at-age (𝑝𝑚𝑎𝑡𝑎 ) was calculated as: ∑𝐼𝑖 𝛾𝑖 ∗ 𝑃(𝐴)365(𝑎−1)+𝑠𝑝,𝑖 ∗ 𝑝𝑚𝑎𝑡365(𝑎−1)+𝑠𝑢,𝑖 (7.6) 𝑝𝑚𝑎𝑡𝑎 = ∑𝐼𝑗 𝛾𝑗 ∗ 𝑃(𝐴)365(𝑎−1)+𝑠𝑝,𝑖 This takes the average proportion mature for the day on which the maturity schedule was calculated (su = 227, or August 15th, also obtained from the modeling subcommittee models) in each year across individual growth morphs, weighted by a growth morph’s probability of being in the population (𝛾𝑖 ) and its probability of being alive on the day of spawning in a year. This calculation produces an average schedule of maturity for each age. 67 Yield-per-recruit analyses To create YPR curves, we calculated YPR for a series of values of F ranging from 0 to 2. For the standard approach, the discrete-time generic method was used (Quinn & Deriso, 1999a), as outlined in Eq. 8.1-8.6. First, 𝐹𝑎 was calculated using the estimated selectivity pattern 𝑠𝑒𝑙𝑎 and baseline fully selected fishing mortality F=0.1 (Eq. 8.1). Then, annual natural mortality was added to 𝐹𝑎 to obtain the total mortality rate 𝑍𝑎 (Eq. 8.2). Initial abundance was set to 1 to make the resulting calculations per-recruit. Abundance at ages above 1, 𝑁𝑎>1, were calculated using the exponential mortality model (Eq. 8.3), and the Baranov catch equation was used to obtain catch-at-age 𝐶𝑎 (Eq. 8.4). This was multiplied by the average weight-at-age during harvest, 𝑊𝑎𝐻 , to obtain yield (Eq. 8.5) and then was summed across all ages to calculate the total YPR (Eq. 8.6). 𝐹𝑎 = 𝐹 ∗ 𝑠𝑒𝑙𝑎 (8.1) 𝑍𝑎 = 𝐹𝑎 + 𝑀 (8.2) 𝑁𝑎+1 = 𝑁𝑎 ∗ 𝑒 −𝑍𝑎 (8.3) 𝐹𝑎 𝐶𝑎 = ∗ 𝑁𝑎 ∗ (1 − 𝑒 −𝑍𝑎 ) 𝑍𝑎 (8.4) 𝑌𝑎 = 𝐶𝑎 ∗ 𝑊𝑎𝐻 (8.5) 50 𝑌𝑃𝑅 = ∑ 𝑌𝑎 𝑎=1 (8.6) True YPR that accounts for growth variability and size selective mortality was calculated using each super-individual’s probability of being alive and probability of harvest using Eq. 6.1-6.4 for a given value of 𝐹. To obtain yield for a super-individual, the probability of harvest on a given day (𝑃(𝐻)𝑑,𝑖 ) was multiplied by its weight on that day (𝑊 ̅𝑑,𝑖 ). These values were then summed for each growth morph across all days (getting the lifetime yield for a super-individual) and 68 weighted as the average of those results using each growth morph’s weight 𝛾. The YPR curves were plotted for comparison. Spawning stock biomass-per-recruit analyses To compare SSBR curves for each method, SSBR was calculated for the same series of F values used in YPR analyses. The standard approach to calculating SSBR was also based on the discrete-time generic method (Quinn and Deriso, 1999). Abundance-at-age was calculated using Eq. 8.1-8.3, and for simplicity, we assumed all individuals were female. First, abundance-at-age on day of spawning, 𝑁𝑠𝑝,𝑎 , was calculated by incorporating the proportion of the year that passed until day of spawning (0.833 is day 304 converted to years, obtained from the 1836 treaty waters Modeling Subcommittee stock assessment models for Lake Superior) (Eq. 9.1). Next, spawning stock biomass, SSBa, was calculated by multiplying 𝑁𝑠𝑝,𝑎 by the average weight-at-age on day of spawning (Eq. 7.5) and proportion mature-at-age on day of spawning (Eq. 7.6) (Eq. 9.2). Finally, SSBa was summed across all ages to obtain SSBR (Eq. 9.3). 𝑁𝑠𝑝,𝑎 = 𝑁𝑎 ∗ 𝑒 − 0.833(𝑍𝑎) (9.1) 𝑆𝑆𝐵𝑎 = 𝑁𝑠𝑝,𝑎 ∗ 𝑊𝑎𝑠𝑝 ∗ 𝑃𝑚𝑎𝑡𝑠𝑝,𝑎 (9.2) 50 𝑆𝑆𝐵𝑅 = ∑ 𝑆𝑆𝐵𝑎 (9.3) 𝑎=1 The true SSBR for a given F value was calculated using each super-individual’s probability of being alive on each day (Eq. 6.1-6.2). From this matrix the probability of being alive on the day of spawning was extracted for each year and individual. For each individual, weight at the beginning of the day of spawning (or 𝑊𝑑=𝑠𝑝,𝑖 ) and the probability of being mature on the day of the survey (August 15th) were extracted from maturity schedule parameters were calculated based on lengths at the time of that survey (Shawn Sitar, Michigan DNR, personal 69 communication). For each individual, the probability of being alive on day of spawning was multiplied by the weight on day of spawning and the probability of being mature, summed across each individual’s potential lifespan (to obtain an individual-specific SSBR value), and weighted using 𝛾. The SSBR curves were plotted for comparison. Reference points Reference points were calculated using standard and true approaches to YPR and SSBR. The reference points included FMAX (the value of F that maximizes YPR), F0.1 (the value of F whose slope produces 10% of the YPR curve when F = 0), and FX% (the value for F that produces X% of the unfished SSBR for X = 35% and X=50%). 70 RESULTS Simulated growth Back-transformed midpoints of VBGM parameter distributions ranged over 480.0-1103.4 for L∞, 0.05-0.17 for K, and 74.6-149.6 for L1 (Figure 10, Table 8). Standardized weights (𝛾) for each growth morph ranged from 3.093e-12 to 4.457e-03} because some combinations of VBGM parameters were far more likely to occur in the Lake Superior system from which we obtained the distributional parameters. Different combinations produced very different growth curves for individuals and the interaction between growth and length-based selective mortality was evident in selectivity patterns. Fish representing growth morphs with low L∞ sometimes failed to reach full selectivity, whereas fish from growth morphs with higher L∞ values reached full selectivity before ~age 14 (Figure 11). Calculations of 𝑁𝑎 , 𝐶𝑎 , 𝑠𝑒𝑙𝑎 , 𝑊𝑎𝐻 , 𝑊𝑎𝑠𝑝 , and 𝑝𝑚𝑎𝑡𝑎 𝑁𝑎 started at 1 and decreased to an asymptote near 0 as age increased due to exponential mortality (Figure 12). By age 30, several population parameters reached an apparent asymptote: 𝑁𝑎 decreased to ~0.0004 by age 50 and 𝐶𝑎 peaked at 0.0115 at a = 13 before. 𝐶𝑎 was dome- shaped and then decreased to near 0 as numbers of fish at higher ages that were large enough to be susceptible to fishing gear became scarcer (Figure 13). 𝐶𝑎 as a proportion of 𝑁𝑎 peaked at 7.6% of 𝑁𝑎 at age 33 (Figure 14). Selectivity based on 𝑁𝑎 and 𝐶𝑎 reached an asymptote near 1 around age 30 (Figure 15). Selectivity-at-age decreased at the highest ages, likely because the method used to calculate selectivity-at-age (Eq. 6) was not able to account for individual-specific selectivity for growth morphs. Average proportion mature-at-age peaked at 0.83 at age 34 (Figure 16). Average weight-at-age increased to a maximum at 2851.4 g at age 44 before 71 decreasing slightly at the highest ages (Figure 17). This decrease at older ages is most likely due to the effects of size-selective removal of large fish earlier in their lifespans, resulting in the loss of fish with larger asymptotic lengths. Similarly, average weight-at-age at time of spawning decreased slightly after a peak of 2579.5 g at age 39 (Figure 17). The calculated average proportion mature-at-age never reached 1 (all fish were not mature at any age). This is in part due to the parameters we used to characterize individual variability in growth and the maturity schedule: one of the 10 L∞ values used for the growth morphs was lower than the inflection point of the maturity schedule (𝐿50 = 558.08 mm). Standard vs. true SSBR and YPR curves The YPR curves calculated with the standard approach predicted lower values than the true approach when F < 0.1, higher values when F > 0.1, and similar values when F = 0.1, where the analyst’s data for 𝐶𝑎 and 𝑊𝑎𝐻 were calculated with a baseline mortality of F = 0.1 (Figure 18). As F increased above 0.1, the curves had similar shapes. At F = 0.5, 1.0, and 2.0, the standard method predicted YPR of 0.062, 0.082, and 0.092 kg higher than true values corresponding to 119%, 126%, and 132% of the true values. The standard SSBR curve predicted lower SSBR for values of F < 0.35 (Figure 19). When F exceeded 0.35, the standard approach slightly overestimated SSBR. The standard method predicted SSBR 0.035 kg higher at F = 0.5 (114%), 0.043 kg higher at F = 1.0 (1185%), and 0.019 kg higher at F = 2.0 than the true method (at F = 2.0, both methods predict numbers very close to 0). The two curves had different intercepts at F = 0 and the standard method predicted SSBR = 3.78 kg lower (67.4%) of the true method intercept = 5.61 kg. 72 Standard method estimates vs. true reference points All YPR-based biological reference points (BRPs) calculated by the standard method were overestimated (Table 9). The standard FMAX value exceeded the true value by 0.498, or 242%. The standard F0.1 estimates were closer, only exceeding the true value by 0.044, or 131%. Each of these YPR-based BRPs had higher YPR predicted for the standard method and exceeded the true YPR by 0.070 kg (122%) for FMAX and 0.039 kg (114%) for F0.1. At values lower than F = 0.35, the standard method underestimated SSBR by not accounting for individuals that matured quickly at younger age or reached greater weight at younger age than others. Biological reference points were 0.026 higher than the true values for F 35% and 0.015 higher for F50% (Figure 11). The standard method underestimated SSBR at low values of F, so reference points F35% and F50% produced SSBR values 67% lower than true SSBR. 73 DISCUSSION This study builds on a body of research investigating individual variability in growth and its effect on estimation in population models and management advice (Kraak et al., 2019a; Kristiansen & Svåsand, 1998; Lowerre-Barbieri et al., 1998; Parma & Deriso, 1990; Ricker, 1969a; Stawitz et al., 2019). Size-selective mortality can preferentially remove faster growing individuals, and if unaccounted for, this could lead to underestimation of spawning stock biomass (Ricker, 1969b). Since then, others have shown that ignoring this variation in models can bias predictions of individual growth trajectories. Sainsbury (1980) found that most population models assume a mean value for size-at-age is representative of an age class, and this remains a common assumption. Models that do not account for individual variation can bias growth parameter estimates across multiple estimation methods (Pilling et al., 2002). Further, these biases can be propagated into stock assessment models or other models used to calculate management reference points (Parma & Deriso, 1990). Our model incorporated individual variability in growth by assuming that individuals varied in all three parameters of the von Bertalanffy growth model. Many studies generally focus on individual variability in L∞ and K, and rarely in the third parameter, which is often t0 (time when the organism length is 0) (Lowerre-Barbieri et al., 1998; Sainsbury, 1980; Shelton & Mangel, 2012) , although Martinez-Garmendia 1998 treated L∞, K, and L1 as varying among individuals (Martínez-Garmendia, 1998). Studies of individual variability have used length-frequency data (Wang & Ellis, 1998), tagging data (Kirkwood and Somers, 1984; Sainsbury, 1980; Wang and Thomas, 1995), and back-calculated increment data (Pilling et al., 2002). Back-calculation and its interaction with size-selective mortality can bias growth models and their interactions with population models (Ricker, 1969b; Schirripa & Trexler, 2000), though some of these effects can be 74 lessened by using a mixed model that treats longitudinal data as a random effect (Escati-Peñaloza et al., 2010; Pilling et al., 2002). Our estimates of growth parameters from a study that treated individual deviations in length increments as random effects assumed that no serious biases were propagated forward from back-calculation (Chapter 1). Our study found bias in per-recruit model predictions when individual variability in growth and its interaction with size-selective mortality was ignored, which is substantial enough to meaningfully impact reference points. FX%-type reference points (where X represents the target or limit percentage of unfished biomass per recruit) have been used to determine whether estimated levels of mortality are too high given the life history of the fish (Gabriel & Mace, 1999; Lenart & Caroffino, 2019). For example, managers using F35% as a limit fishing intensity are assuming that this reference point corresponds to a value of spawning stock biomass they wish to preserve (35% of the unfished SSBR). In our study, the standard method produced a value of F 25% greater than the true value; therefore, a manager using reference points calculated via the standard method has a biased perception of risk to the population and in the case outlined here would be recommending too high a fishing intensity for their target SSBR. Variability in growth, when ignored, causes overestimation of reference points like F 0.1 and overestimation of SSBR and YPR (especially at higher levels of F), and that incorporating stochasticity in individual growth parameters reduces this bias (Lowerre-Barbieri et al., 1998; Parma & Deriso, 1990). We found that a standard method that ignored individual variability underestimated SSBR at low levels of F. This could be attributed to some individuals being larger and more likely to be mature faster than would be predicted with mean length-at-age models, thereby contributing to spawning stock by weight earlier than mean weight at age would predict. Further, average weight-at-spawning calculated in Eq. 7.5, which is used in the standard 75 method, results in lower weights across all ages than is actually the case under baseline mortality (Figure 20). The average proportion mature-at-age used in the standard method shows a dip at higher ages, which in the standard SSBR method resulted in fewer spawners with slightly lower weights, even at F = 0, because the standard method curves are based on “observed values” for selectivity and maturity for which there was an effect of fishing. Our model found that YPR and SSBR were overestimated when individual variability was ignored at higher levels of F. We modeled individual variability as consistent over the lifespan of individuals with no time- varying parameters. Key differences between our model and the Parma & Deriso (1990) framework, which uses age and size distributions rather than individual growth trajectories, are in the sources of variation: their variability in growth was attributed entirely to environmental effects and allowed for parameters to change over time due to size-selective mortality. Further, they assumed that fishing mortality and natural mortality occurred at separate points during the year. Seasonal growth variability and individual produced higher YPR and egg production than without seasonal growth variability (Lowerre-Barbieri et al., 1998). Our model does not assume seasonality in growth, although lake trout growth, maturation, and habitat vary seasonally (Cline et al., 2013; Goetz et al., 2011; Kao et al., 2015; Lawrie, 1963; von Biela et al., 2021). Our growth model could be modified to incorporate seasonal variability (i.e., using an oscillating von Bertalanffy model), but was outside the scope of our analysis. However, several studies investigating individual variability in growth in per-recruit models have shown that temporal variability of growth parameters contributes to bias (Miller et al., 2018; Stawitz et al., 2019; You- Gan Wang & Thomas, 1995). The relative contribution of intrinsic growth variation (wherein the individual-specific parameters are constant across their lifespan) and extrinsic sources of variation (i.e., annual variability in temperature, environmental effects, etc.), also referred to as 76 persistent and transient variation [sensu (Webber & Thorson, 2016)], can be extremely difficult to estimate from standard fishery data (Parma & Deriso, 1990). The model from which we obtained our growth parameter distributions, and the growth model we used, incorporated both persistent and transient variability in the form of individual-specific growth parameters and a year- and individual-specific process error (Chapter 1). These parameter estimates suggested that for our study system, transient variability did not contribute as much to variation in length-at-age as persistent variation. The present study shows that persistent variation is substantial and, when ignored, can result in biased estimates of biological reference points and per-recruit model output. We focused on a basic formulation of a per-recruit model (Gabriel et al., 1989) and used life history parameters from Great Lakes lake trout statistical catch-at-age (SCAA) models (Caroffino & Lenart, 2000). Lake trout harvest policies in the Great Lakes incorporate output from such SCAA models, which currently use mean weight-at-age in the calculation of spawning stock (Shawn Sitar, Michigan DNR, personal communication). Though we do not perform a full assessment, our model simulated a Lake Superior population with growth variation based on estimates from that system and calculated reference points. Therefore, our results should be of interest to harvest policies or population models that are modeling, in any capacity, growth of lake trout or using metrics of growth like length- or weight-at-age. Our model, which incorporates individual variability in growth and its interaction with size-selective mortality, to calculate probabilities of survival and harvest, shows that YPR, SSBR, and biological reference points can be mis-estimated when individual variability is ignored. We recommend that managers who use these models in the Great Lakes, or similar models for different systems, consider either accounting for individual variability in growth (using our proposed model or one 77 of several other methods outlined above), or acknowledge the potential bias in reference points and prescribe conservative policy that accounts for the bias. Our methods have several limitations. To examine specifically the influence of individual variability in growth, we restricted uncertainty in other areas of the model and made several limiting assumptions. Firstly, we did not account for seasonality in lake trout growth and instead treated within-year growth as constant, although the influence of seasonality in growth for a slow-growing and long-lived species such as the lake trout may be less important. We allowed for year-to-year fluctuations when calculating growth increments via a random process error, but the standard deviation of that process error was small compared to VBGM parameter standard deviations. Secondly, we fixed background natural mortality at 0.1 as constant over time. This level of natural mortality seems reasonable for managed lake trout units in the Great Lakes (Lenart & Caroffino, 2019) but still carries implicit assumptions about conditions in Lake Superior for lake trout. While we assessed the impact of fishing mortality via selectivity curves and individual variation, we did not explicitly model mortality by the invasive Sea Lamprey (Petromyzon marinus), which has been found to be size-selective in predation (Bence et al., 2003; Sitar & He, 2006). Given that sea lamprey attacks are generally believed to follow a logistic function and at larger sizes(Bence et al., 2003), it is plausible that the results from this study would apply to the consequences of changes in sea lamprey mortality even if fishing mortality stayed at a low level. We also fixed the maturation schedule according to existing parameterizations used by management population models for Lake Superior lake trout, to limit uncertainty in other areas of the model. Variability in these areas could be incorporated into our framework, but the aim of this study was to demonstrate the bias that could exist in a simplified system with individual variability in growth and we were concerned about compromising our 78 ability to cleanly interpret results. The patterns modeled here, including growth function, maturation and selectivity patterns, and weight-length relationships, can be easily modified for different life histories and exploitation patterns. This said, we have no reason to believe that various complexities and increased realism would alter the fundamental result, that given the nature of among individual variation in growth in lean lake trout in Lake Superior, ignoring this variation has a real potential to bias management reference points and put fishery management objectives at risk. Our analysis also assumed that the analyst had access to accurate data for abundance-at-age, catch-at-age, and proportion mature-at-age, as a best-case scenario for the standard method. Some studies have found that bias caused by uncertainty in bad data can exceed bias caused by individual variability (Stawitz et al., 2019). We also simulated in daily increments, which offers a specificity in incremental data that may not be possible for temperate species like lake trout but could be helpful if applying this model to fish with shorter lifespans, like in tropical regions. Fourthly, we assume a simple logistic selectivity pattern. The SCAA model used for lake trout in Lake Superior assumes a double logistic for gill nets, and thus depending on the type of fishery, dome-shaped selectivity may also be appropriate (Hansen et al., 1997). Our model framework could be modified to incorporate a more complex selectivity pattern, and this could have significant effects on results. Fifthly, values of L∞ we used meant that not all individuals reached maturity in their lifespans because of our choice of L50, and because this parameter was not linked to growth trajectories of individual growth nodes. This is highly unrealistic and must necessarily be modified before application of the model to management. One possible approach would be to assume that L50 and A50 lie on a “reactive norm” curve (Breck et al., 2020). Past work for other species has suggested that slower growing fish mature at an older age but smaller size (Angilletta et al., 2004; Heino et al., 2002; Hutchings & 79 Jones, 1998; Shuter et al., 2016). The reactive norm idea might feasibly allow the calculation of morph specific inflection points for the logistic maturity function. We also assumed that the standard approach to calculating spawning biomass per recruit would use weight at age based on all fish in the population. Our experience is that this is what is often done, but an analyst could calculate weight at age only for mature fish for use in calculations of spawning biomass, which would remove some of the bias at low levels of F for spawning biomass per recruit estimates. Finally, we fixed the exponent in our weight-length relationship at 3.19 and the coefficient at 2.35e-06, both estimated from our dataset. However, the value of this exponent varies for lake trout across different regions in the Great Lakes (He et al., 2020; He & Bence, 2007). Changing the exponent and coefficient can significantly impact YPR and SSBR estimates, especially at higher lengths, and further work would need to assess sensitivity of our results to this variable. Our study is based on a model using one specific morphotype of lake trout (lean). Lake trout morphs can vary in habitat use and shape (Muir et al., 2014), as well as life history parameters (Hansen, Nate, Muir, et al., 2016), so our results are only applicable to the lean morphotype. Despite these limitations, results showed that populations where significant variation in growth parameters has been measured can produce biased reference points and per-recruit calculations if this variability is ignored. Our analyses made assumptions about where variability is included, but showed overestimation of reference points, yield, and biomass, especially at higher levels of fishing intensity. Because our model simulates in daily time steps, it can be modified for fish with very short lifespans. Further, the choice of growth model can be easily changed depending on the focus of the study and other components (i.e., selectivity, maturity) can be modified to reflect conditions of the system. Ultimately, we showed that with lake trout in the Great Lakes, 80 models used for management may overestimate the potential of the stock and we recommend investigating the feasibility of incorporating individual variability into existing models. 81 APPENDIX 82 APPENDIX Table 7 Name Symbol Description Units Assigning VBGM parameter values L infinity L∞ VBGM parameter. Asymptotic length mm VBGM parameter. Brody growth K K coefficient. 1/years L1 L1 VBGM parameter. Length at age 1 mm Midpoint of each bin (10 total) within a Midpoint mp VBGM parameter distribution mm Lower limit calculated for each bin (10 total) within a VBGM parameter Bin lower limit a distribution mm Upper limit calculated for each bin (10 total) within a VBGM parameter Bin upper limit b distribution mm Simulating growth Individual- and age-specific process Process error δ error Process error standard Standard deviation of random variable deviation σ_ δ δ, process error Increment ΔL Increment in length mm Weight-length relationship Body weight w Body weight of a fish kg Body length L Body length of a fish mm Coefficient of weight-length Weight-length coefficient λ relationship Weight-length exponent β Exponent of weight-length relationship Growth morph weighting Probability density function value for each growth morph's specific combination of VBGM, calculated using the estimated trivariate distribution from Weighting factor γ the previous chapter Table 7: Parameter and variable names, symbols, descriptions, and units (where appropriate). 83 Table 7 (cont’d) Selectivity Pattern Selectivity-at-length, day sel Selectivity-at-length, at a given time Length of a fish on a given day, used to Length L calculate its selectivity for that day mm Length of a fish at entry to the fishery Length-at-entry LE (inflection point in the logistic curve) mm Alpha α Slope of the logistic function Maturity pattern Probability of being mature- Probability of being mature at-length, at at-length, day pmat_(L,d) a given time Length of a fish on a given day, used to calculate its probability of being mature Length-at-day Ld for that day mm Length of a fish at which 50% are mature, or the inflection point for the L50 L50 logistic function mm Beta β Slope of the logistic function Population dynamics Instantaneous fishing mortality F Instantaneous fishing mortality Natural mortality M Natural mortality Total mortality Z Total mortality Abundance/Probability of being alive N or P(A) Abundance, or probability of being alive Probability of harvest P(H) Probability of being harvested Catch C Catch (in numbers) Average weight at midpoint of day 𝑤 ̅ Average weight at midpoint of day kg Growth G Geometric growth Weight-at-age at time of Average weight-at-age at time of harvest Wa,H harvest kg Weight-at-age at time of Average weight-at-age at time of spawning Wa,sp spawning kg Per-recruit models Yield Y Yield (in weight) kg Yield (in weight) produced over the life Yield-per-recruit YPR of an individual fish kg Spawning stock biomass SSB Spawning stock biomass (in weight) kg Spawning stock biomass (in weight) Spawning stock biomass- produced over the life of an individual per-recruit SSBR fish kg 84 Table 7 (cont’d) Subscripts Day d Day of year Age, or year (spanning 1:50 for this Age/year a model) Sometimes j is used in weighted average Individual i calculations Day of spawning sp Day on which spawning occurs Day on which the summer survey occurs, creating parameters for maturity Day of summer survey su schedules Harvest H Denotes occurring on the day of harvest Table 8 L infinity K L1 Log Log Log scale Backtransformed scale Backtransformed scale Backtransformed 6.17 480.02 -2.92 0.05 4.31 74.59 6.34 566.77 -2.69 0.07 4.45 85.71 6.43 618.86 -2.56 0.08 4.52 92.24 6.50 663.47 -2.47 0.08 4.58 97.77 6.56 706.16 -2.38 0.09 4.63 103.01 6.62 750.06 -2.30 0.10 4.69 108.33 6.68 798.32 -2.21 0.11 4.74 114.13 6.75 855.87 -2.12 0.12 4.80 120.97 6.84 934.53 -1.99 0.14 4.87 130.20 7.01 1103.42 -1.76 0.17 5.01 149.60 Table 8: von Bertalanffy parameter values used. 10 values of each von Bertalanffy parameter representing a range of potential values as estimated for Lake Superior lake trout. Values are presented on log scale and back-transformed. The 10 values of each parameter were combined into 1000 unique triplet nodes which we refer to as “growth morphs”. 85 Table 9 Biological Relative Reference True Standard Relative True Standard YPR/SSBR Point F F F Bias YPR/SSBR YPR/SSBR bias FMAX 0.349 0.847 -1.425 0.321 0.392 -0.219 F0.1 0.142 0.187 -0.312 0.285 0.324 -0.138 F35% 0.107 0.133 -0.248 1.963 1.323 0.326 F50% 0.064 0.080 -0.241 2.805 1.890 0.326 Table 9: Biological reference points. Biological reference points (BRPs) calculated using standard and true methods. The value of F and the associated yield-per-recruit (YPR) or spawning stock biomass-per-recruit (SSBR) values are shown. Relative bias was calculated by subtracting the standard estimate from the true estimate and dividing by the true estimate. 86 Figure 9: Individual length-at-age of Lean Lake Trout. Individual length-at-age of lake trout back-calculated from ototlith data collected from n=410 individuals in Lake Superior. These data were used to estimate individual variability in growth. 87 Figure 10: Growth parameter values used to simulate growth morphs. Distribution of L∞, K, and L1 values using the estimated log-scale mean and numerical-scale standard deviation estimated in the previous chapter. Midpoints are shown in red dashed lines. 88 Figure 11: Individual-specific selectivity patterns and individual-specific growth patterns. Examples of super-individual-specific selectivity patterns and growth curves calculated using the individual’s length-at-day and the respective von Bertalanffy growth parameters. Standardized weights (probability of observing that growth morph’s von Bertalanffy parameters) of each growth morph (from top of legend to bottom) are 0.00445, 0.00053, 0.00076, and 0.00045. 89 Figure 12: Numbers-at-age. Numbers-at-age begin at 1, meaning values represent proportions, or per-recruit values. Figure 13: Catch-at-age. Catch-at-age per recruit calculated from simulated data. 90 Figure 14: Catch-at-age as a proportion of numbers-at-age (i.e., exploitation rate). Figure 15: Selectivity-at-age. Selectivity-at-age pattern calculated using catch-at-age per recruit and numbers-at-age per recruit assuming the analyst has good data. 91 Figure 16: Proportion mature-at-age. Average proportion mature-at-age calculated from simulated data. The curve never reaches 1 because some fish never reach the size of maturity due to the variability we used for L∞. Figure 17: Weight-at-age (spawning & harvest). Average weight-at-age at time of harvest (black) and at time of spawning (red) calculated from simulated data. 92 Figure 18: Estimated YPR curves. YPR curves, calculated using the standard method compared to true values. Vertical lines represent biological reference points F MAX and F0.1 calculated for the standard and true methods. Figure 19: Estimated SSBR curves. Standard method estimates and true SSBR curves with biological reference points with increased resolution. Vertical lines represent biological reference points F35% and F50% for both standard and true methods. The curves stay converged and get increasingly close at F > 0.3 to F = 2. 93 Figure 20: Average weight-at-age at different levels of F. Average weight-at-age at harvest, calculated as was done for the standard method for different levels of fishing intensity F. Colors represent different values of F. The standard methods use average weight-at-age calculated at a baseline level of mortality (F = 0.1), though as F changes the average weight-at-age declines. 94 LITERATURE CITED 95 LITERATURE CITED Angilletta, M.J., Steury, T.D., Sears, M.W., 2004. Temperature, Growth Rate, and Body Size in Ectotherms: Fitting Pieces of a Life-History Puzzle 1. Bence, J.R., Bergstedt, R.A., Christie, G.C., Cochran, P.A., Ebener, M.P., Koonce, J.F., Rutter, M.A., Swink, W.D., 2003. Sea lamprey (Petromyzon marinus) parasite-host interactions in the Great Lakes. Journal of Great Lakes Research 29, 253–282. https://doi.org/10.1016/S0380-1330(03)70493-6 Breck, J.E., Simon, C.P., Rutherford, E.S., Low, B.S., Lamberson, P.J., Rogers, M.W., 2020. The geometry of reaction norms yields insights on classical fitness functions for Great Lakes salmon. PLoS ONE 15. https://doi.org/10.1371/journal.pone.0228990 Burkardt, J., 2014. The Truncated Normal Distribution [WWW Document]. Department of Scientific Computing, Florida State University. URL http://people.sc.fsu.edu/∼jburkardt/presentations/truncated normal.pdf (accessed 6.27.22). Caroffino, D.C., Lenart, S.J., 2000. Statistical catch-at-age models used to describe the status of lean lake trout populations in the 1836-Treaty ceded waters of lakes Michigan , Huron , and Superior at the inception of the 2000 Consent Decree A Report Completed by the Modeling Subcommittee. Cline, T.J., Bennington, V., Kitchell, J.F., 2013. Climate Change Expands the Spatial Extent and Duration of Preferred Thermal Habitat for Lake Superior Fishes. PLoS ONE 8. https://doi.org/10.1371/journal.pone.0062279 Escati-Peñaloza, G., Parma, A.M., Orensanz, J.M.L., 2010. Analysis of longitudinal growth increment data using mixed-effects models: Individual and spatial variability in a clam. Fisheries Research 105, 91–101. https://doi.org/10.1016/j.fishres.2010.03.007 Fabens, A.J., 1965. Properties and fitting of the von Bertalanffy Growth Curve. Growth 29, 265– 289. Flinn, S.A., Midway, S.R., 2021. Trends in growth modeling in fisheries science. Fishes 6, 1–18. https://doi.org/10.3390/fishes6010001 Gabriel, W.L., Mace, P.M., 1999. A review of biological reference points in the context of the precautionary approach. Proceedings of the fifth national NMFS stock assessment workshop: providing scientific advice to implement the precautionary approach under the Magnuson- Stevens fishery conservation and management act. NOAA Tech Memo NMFS-F/SPO-40 34–45. Gabriel, W.L., Sissenwine, M.P., Overholtz, W.J., 1989. Analysis of Spawning Stock Biomass per Recruit: An Example for Georges Bank Haddock. North American Journal of Fisheries 96 Management 9, 383–391. https://doi.org/10.1577/1548- 8675(1989)009<0383:aossbp>2.3.co;2 Genz, A., Bretz, F., Miwa, T., Mi, X., Leisch, F., Scheipl, F., Hothorn, T., 2020. mvtnorm: Multivariate Normal and t Distributions. Goetz, F., Sitar, S., Rosauer, D., Swanson, P., Bronte, C.R., Dickey, J., Simchick, C., 2011. The reproductive biology of Siscowet and lean lake trout in southern lake superior. Trans Am Fish Soc 140, 1472–1491. https://doi.org/10.1080/00028487.2011.630276 Hansen, M.J., Madenjian, C.P., Selgeby, J.H., Helser, T.E., 1997. Gillnet selectivity for lake trout (Salvelinus namaycush ) in Lake Superior. Canadian Journal of Fisheries and Aquatic Sciences 54, 2483–2490. https://doi.org/10.1139/f97-156 Hansen, M.J., Nate, N.A., Muir, A.M., Bronte, C.R., Zimmerman, M.S., Krueger, C.C., 2016. Life history variation among four lake trout morphs at Isle Royale, Lake Superior. Journal of Great Lakes Research 42, 421–432. https://doi.org/10.1016/j.jglr.2015.12.011 He, J.X., Bence, J.R., 2007. Modeling Annual Growth Variation using a Hierarchical Bayesian Approach and the von Bertalanffy Growth Function, with Application to Lake Trout in Southern Lake Huron. Trans Am Fish Soc 136, 318–330. https://doi.org/10.1577/t06-108.1 He, J.X., Bence, J.R., Johnson, J.E., Clapp, D.F., Ebener, M.P., 2008. Modeling Variation in Mass-Length Relations and Condition Indices of Lake Trout and Chinook Salmon in Lake Huron: A Hierarchical Bayesian Approach. Trans Am Fish Soc 137, 801–817. https://doi.org/10.1577/t07-012.1 He, J.X., Bence, J.R., Madenjian, C.P., Claramunt, R.M., 2020. Dynamics of lake trout production in the main basin of Lake Huron. ICES Journal of Marine Science 77, 975–987. https://doi.org/10.1093/icesjms/fsaa030 Heino, M., Dieckmann, U., Godø, O.R., 2002. Measuring Probabilistic Reaction Norms for Age and Size at Maturation, Evolution. Hutchings, J.A., Jones, M.E.B., 1998. Life history variation and growth rate thresholds for maturity in Atlantic salmon, Salmo salar. Canadian Journal of Fisheries and Aquatic Sciences 55, 22–47. Kao, Y.C., Madenjian, C.P., Bunnell, D.B., Lofgren, B.M., Perroud, M., 2015. Temperature effects induced by climate change on the growth and consumption by salmonines in Lakes Michigan and Huron. Environmental Biology of Fishes 98, 1089–1104. https://doi.org/10.1007/s10641-014-0352-6 Kirkwood, G.P., Somers, I.F., 1984. Growth of two species of tiger prawn, penaeus esculentus and p. Semisulcatus, in the western gulf of carpentaria. Marine and Freshwater Research 35, 703–712. https://doi.org/10.1071/MF9840703 97 Kraak, S.B.M., Haase, S., Minto, C., Santos, J., Anderson, E., 2019a. The Rosa Lee phenomenon and its consequences for fisheries advice on changes in fishing mortality or gear selectivity. ICES Journal of Marine Science 76, 2179–2192. https://doi.org/10.1093/icesjms/fsz107 Kraak, S.B.M., Haase, S., Minto, C., Santos, J., Anderson, E., 2019b. The Rosa Lee phenomenon and its consequences for fisheries advice on changes in fishing mortality or gear selectivity. ICES Journal of Marine Science 76, 2179–2192. https://doi.org/10.1093/icesjms/fsz107 Kristiansen, T.S., Svåsand, T., 1998. Effect of size-selective mortality on growth of coastal cod illustrated by tagging data and an individual-based growth and mortality model. Journal of Fish Biology 52, 688–705. https://doi.org/10.1006/jfbi.1997.0611 Kuparinen, A., Kuikka, S., Merilä, J., 2009. Estimating fisheries-induced selection: Traditional gear selectivity research meets fisheries-induced evolution. Evolutionary Applications 2, 234–243. https://doi.org/10.1111/j.1752-4571.2009.00070.x Lawrie, A.H., 1963. Intraseasonal Growth in Lake Superior Lake Trout. Journal of the Fisheries Research Board of Canada 20, 491–496. https://doi.org/10.1139/f63-037 Lee, R.M., 1912. An investigation into the methods of growth determination in fishes by means of scales. ICES Journal of Marine Science s1, 35. Lenart, S.J., Caroffino, D.C., 2019. Technical Fisheries Committee Administrative Report 2019: Status of Lake Trout and Lake Whitefish Populations in the 1836 Treaty-Ceded Waters of Lakes Superior, Huron and Michigan, with Recommended Yield and Effort Levels for 2019. Lowerre-Barbieri, S.K., Lowerre, J.M., Barbieri, L.R., 1998. Multiple spawning and the dynamics of fish populations: Inferences from an individual-based simulation model. Canadian Journal of Fisheries and Aquatic Sciences 55, 2244–2254. https://doi.org/10.1139/f98-105 Martínez-Garmendia, J., 1998. Simulation analysis of evolutionary response of fish populations to size-selective harvesting with the use of an individual-based model. Ecological Modelling 111, 37–60. https://doi.org/10.1016/S0304-3800(98)00093-3 Miller, T.J., Brien, L.O., Fratantoni, P.S., 2018. Temporal and environmental variation in growth and maturity and effects on management reference points of Georges Bank. Canada Journal of Fisheries and Aquatic Sciences 75, 2159–2171. Muir, A.M., Bronte, C.R., Zimmerman, M.S., Quinlan, H.R., Glase, J.D., Krueger, C.C., 2014. Ecomorphological Diversity of Lake Trout at Isle Royale, Lake Superior. Trans Am Fish Soc 143, 972–987. https://doi.org/10.1080/00028487.2014.900823 Parma, A.M., Deriso, R.B., 1990. Dynamics of Age and Size Composition in a Population Subject to Size-Selective Mortality: Effects of Phenotypic Variability in Growth. Canada Journal of Fisheries and Aquatic Sciences 47. 98 Pilling, G.M., Kirkwood, G.P., Walker, S.G., 2002. An improved method for estimating individual growth variability in fish, and the correlation between von Bertalanffy growth parameters. Canadian Journal of Fisheries and Aquatic Sciences 59, 424–432. https://doi.org/10.1139/f02-022 Quinn, T.J., Deriso, R.B., 1999a. Age-structured Models: Per-Recruit and Year-Class Models, in: Quantitative Fish Dynamics. Oxford University Press, New York, pp. 239–267. Quinn, T.J., Deriso, R.B., 1999b. Quantitative Fish Dynamics. Oxford University Press, New York. R Core Team, 2020. R: A language and environment for statistical computing. Ricker, W.E., 1969a. Effects of Size-Selective Mortality and Sampling Bias on Estimates of Growth, Mortality, Production, and Yield. Journal of the Fisheries Research Board of Canada 26, 479–541. Ricker, W.E., 1969b. Effects of Size-Selective Mortality and Sampling Bias on Estimates of Growth, Mortality, Production, and Yield. Journal of the Fisheries Research Board of Canada 26. Sainsbury, K.J., 1980. Effect of Individual Variability on the von Bertalanffy Growth Equation. Canadian Journal of Fisheries and Aquatic Sciences 37, 241–247. Schirripa, M.J., Trexler, J.C., 2000. Effects of mortality and gear selectivity on the fish otolith radius-total length relation. Fisheries Research 46, 83–89. https://doi.org/10.1016/S0165- 7836(00)00135-1 Shelton, A.O., Mangel, M., 2012. Estimating von Bertalanffy parameters with individual and environmental variations in growth. Journal of Biological Dynamics 6, 3–30. https://doi.org/10.1080/17513758.2012.697195 Shuter, B.J., Giacomini, H.C., de Kerckhove, D., Vascotto, K., 2016. Fish life history dynamics: Shifts in prey size structure evoke shifts in predator maturation traits. Canadian Journal of Fisheries and Aquatic Sciences 73, 693–708. https://doi.org/10.1139/cjfas-2015-0190 Sitar, S.P., He, J.X., 2006. Growth and Maturity of Hatchery and Wild Lean Lake Trout during Population Recovery in Michigan Waters of Lake Superior. Trans Am Fish Soc 135, 915– 923. https://doi.org/10.1577/t05-019.1 Stawitz, C.C., Haltuch, M.A., Johnson, K.F., 2019. How does growth misspecification affect management advice derived from an integrated fisheries stock assessment model? Fisheries Research 213, 12–21. https://doi.org/10.1016/j.fishres.2019.01.004 99 von Biela, V.R., Black, B.A., Young, D.B., van der Sleen, P., Bartz, K.K., Zimmerman, C.E., 2021. Lake trout growth is sensitive to spring temperature in southwest Alaska lakes. Ecology of Freshwater Fish 30, 88–99. https://doi.org/10.1111/eff.12566 Wang, Y.-G., Ellis, N., 1998. Effect of individual variability on estimation of population parameters from length-frequency data. Canadian Journal of Fisheries and Aquatic Sciences 55, 2393–2401. https://doi.org/10.1139/f98-134 Webber, D.N., Thorson, J.T., 2016. Variation in growth among individuals and over time: A case study and simulation experiment involving tagged Antarctic toothfish. Fisheries Research 180, 67–76. https://doi.org/10.1016/j.fishres.2015.08.016 Wolf, M., Weissing, F.J., 2012. Animal personalities: Consequences for ecology and evolution. Trends in Ecology and Evolution. https://doi.org/10.1016/j.tree.2012.05.001 You-Gan Wang, Thomas, M.R., 1995. Accounting for individual variability in the von Bertalanffy growth model. Canadian Journal of Fisheries and Aquatic Sciences 52, 1368– 1375. https://doi.org/10.1139/f95-132 100