WALLEYE (SANDER VITREUS) DYNAMICS IN THE INLAND WATERWAY, MICHIGAN By Seth J. Herbst A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Fisheries and Wildlife – Doctor of Philosophy 2015 ABSTRACT WALLEYE (SANDER VITREUS) DYNAMICS IN THE INLAND WATERWAY, MICHIGAN By Seth J. Herbst The ecology and population dynamics of fish that inhabit interconnected waterbodies are challenging to understand. Movement of fish among waterbodies presents difficulties in estimating population size and exploitation rate, and movement of native and non-native forage complicates food web dynamics. The Inland Waterway in northern Michigan is a lake chain that consists of four lakes (Burt, Crooked, Mullett, and Pickerel) and connecting rivers that contain walleye, which is a highly mobile species. Walleye within the Inland Waterway are harvested from two distinct fisheries; spearing and recreational-angling fisheries. Harvest quotas are based on spring mark-recapture population estimates of adult walleye for individual lakes using a closed population model and a constant maximum exploitation rate (u = 0.35) harvest control rule. I estimated the post-spawn movement and fishing mortality rates to determine the appropriateness of the current model assumptions for making harvest management decisions. I developed a tag-recovery state-space model fitted using Bayesian estimation techniques to determine the spatial structure of walleye populations in this waterway. Walleye moved extensively among waterbodies, with rates of moving out of a lake or river post spawn ranging between 0.10 and 0.82. This level of intermixing is substantial enough to warrant consideration in the population estimation and harvest allocation processes. I developed a stochastic simulation model to determine the influence of movement patterns on the application of the current harvest control rule. The total allowable catch needed to maintain the desired level of maximum exploitation (u = 0.35) on the individual spawning populations required lake-specific adjustments during the angling harvest period. Although adjustments were needed to meet target exploitation rates, current estimates of exploitation rate were low enough that no individual spawning population was being overexploited. Movement of non-native species into a lake chain system can negatively influence the foraging ecology of native predators. The Inland Waterway has experienced the secondary spread of round goby, alewife, and rainbow smelt from the Great Lakes. Prior to this study, little was known on how these non-native species have been integrated into the foraging pattern of the system’s native top predator. Using stomach contents and stable isotope analysis, I determined that walleye foraging did not always follow the pattern that would be expected based on available literature. We also determined that the non-native round goby became a dominant portion of the walleye diets in Burt and Mullett lakes, which were the two lakes where round goby were abundant. In addition, we determined that walleye did not prey on non-native pelagic fish prey, despite the large area of pelagic habitat and established populations of alewife in Burt and Mullett lakes. This study illustrated that native predators in smaller inland lakes exhibit flexibility in their response to non-native species, but much like predators within the Great Lakes, predators from our smaller scale study lakes have integrated round goby into their forage ecology. ACKNOWLEDGMENTS Special thanks go to the Michigan Department of Natural Resources and the Robert C. Ball and Betty A. Ball Michigan State University Fisheries and Wildlife Fellowship who provided the funding for this project. I would like to thank my graduate advisor, Dr. Daniel Hayes immensely for his contributions to this project, and helping me grow as a professional. Dr. Hayes’ availability, ambition, and intellect have been very influential to my success at Michigan State University and now as a professional with the Michigan Department of Natural Resources (MDNR). His mentoring on research design, data analysis, communication of results, and writing has taught me how to become a better scientist. Dr. Hayes also provided me with a release from the “grind” through introducing me to, and joining me on hunting and fishing trips in the area. I would like to thank my friend and fellow Ph.D. candidate Bryan Stevens. He assisted greatly with model development and ensured that portions of this project were properly coded for analysis. His assistance proved critical during my graduate work by increasing my knowledge of available statistical options for data analysis. I would also like to thank the members of my graduate committee, Dr. Brian Roth, Dr. Mary Bremigan, Dr. Gary Mittelbach, and Nick Popoff for all their valuable input and time on this project, especially with study design, data analysis, and thought process when drafting this document. I would like to thank Patrick Hanchin (MDNR), Tim Cwalinski (MDNR), Neal Godby (MDNR), Pat VanDaele (MDNR), and Maxwell Field (LTTB) for their assistance and willingness to share insight on my study area and also for the long hours in inclement weather that they put in to ensure the data collection process of this project was completed successfully. I iv would like to thank Dr. Jason Stockwell (UVM) and Dr. Jim Bence (MSU) for their assistance with reviewing portions of this document. I would like to thank Ryan MacWilliams, Dan Quinn, Joe Parzych, Mike Rucinski, Elle Gulotty, Kevin Osantowski, and Greg Byford for their hard work and dedication in the field and lab during this project. I am thankful for the rest of my fellow graduate students that assisted with field work. Thanks also to my family and friends for their support during my time as a graduate student at MSU. Who would have thought that fishing trips with my grandparents and parents as a child would spark the passion to earn advanced degrees in fisheries and lead to my current profession?! Most importantly, I am very grateful for my fiancée Liz, for her support and ability to provide much needed distractions from my graduate work during this journey. I am much appreciative of her wedding planning efforts during this chaotic time for me. v TABLE OF CONTENTS LIST OF TABLES ....................................................................................................................... viii LIST OF FIGURES ....................................................................................................................... xi CHAPTER 1 ................................................................................................................................... 1 A State-Space Model for Estimating Walleye (Sander vitreus) Movement and Fishing Mortality in Interconnected Lakes .................................................................................................................. 1 Introduction ................................................................................................................................. 1 Methods ....................................................................................................................................... 5 Study area ................................................................................................................................ 5 Tagging and recovery data ...................................................................................................... 6 Model Structure ....................................................................................................................... 7 General Approach ................................................................................................................... 7 Population process model ....................................................................................................... 8 Tag-recovery observation model ........................................................................................... 11 Prior distributions for model parameters .............................................................................. 12 Model set................................................................................................................................ 15 Model fitting and evaluation.................................................................................................. 16 Results ....................................................................................................................................... 17 Model Selection ..................................................................................................................... 17 Demographic parameters ...................................................................................................... 18 Movement rates...................................................................................................................... 18 Fishing mortality ................................................................................................................... 18 Sensitivity to process error and tag shedding rate assumptions ........................................... 19 Discussion ................................................................................................................................. 21 Walleye Dynamics in the Inland Waterway ........................................................................... 24 CHAPTER 2 ................................................................................................................................. 26 The Influence of Movement Dynamics on Harvest Allocation for Walleye (Sander vitreus) in Intermixed Fisheries in a Chain of Lakes ..................................................................................... 26 Introduction ............................................................................................................................... 26 Methods ..................................................................................................................................... 30 Study area .............................................................................................................................. 30 General approach and model structure................................................................................. 31 Performance metrics.............................................................................................................. 35 Results ....................................................................................................................................... 37 Discussion ................................................................................................................................. 39 Implications of spatial structure on fishery implementation ................................................. 39 Addressing sources of uncertainty ......................................................................................... 42 Management implications and conclusion ............................................................................ 44 vi CHAPTER 3 ................................................................................................................................. 47 Walleye Foraging Ecology in an Interconnected Chain of Lakes Influenced by Non-Native Species .......................................................................................................................................... 47 Introduction ............................................................................................................................... 47 Methods ..................................................................................................................................... 51 Study area .............................................................................................................................. 51 Forage assessment ................................................................................................................. 52 Stomach content analyses ...................................................................................................... 53 Stable Isotope analyses .......................................................................................................... 54 Results ....................................................................................................................................... 57 Forage assessment ................................................................................................................. 57 Stomach content analysis....................................................................................................... 58 Stable isotopes ....................................................................................................................... 59 Discussion ................................................................................................................................. 61 Prey community ..................................................................................................................... 61 Foraging Ecology .................................................................................................................. 62 Stomach content analysis....................................................................................................... 62 Consistency between stomach content and stable isotope analysis ...................................... 63 Relationship between foraging strategy and habitat ............................................................. 64 Management implications of non-natives in the food web .................................................... 65 Conclusion ............................................................................................................................. 67 APPENDICES .............................................................................................................................. 69 APPENDIX 1.0: Supporting tables and figures for Chapter 1: A State-Space Model for Estimating Walleye (Sander vitreus) Movement and Fishing Mortality in Interconnected Lakes. ....................................................................................................................................................... 70 APPENDIX 1.1: The explanation for the derivation of prior distribution for fishing mortality .. 78 APPENDIX 1.2: Number of individuals released and recovered by location and year for each cohort ............................................................................................................................................ 79 APPENDIX 2.0: Supporting tables and figures for Chapter 2: The Influence of Movement Dynamics on Harvest Allocation for Walleye (Sander vitreus) in Intermixed Fisheries in a Chain of Lakes ......................................................................................................................................... 81 APPENDIX 3.0: Supporting tables and figures for Chapter 3: Walleye Foraging Ecology in an Interconnected Chain of Lakes Influenced by Non-Native Species ............................................. 89 APPENDIX 3.1: Stable Isotope summary statistics (δ13C, δ15N, trophic position (TP), and percent littoral). ............................................................................................................................. 96 APPENDIX 3.2: Analysis of covariance results to determine the influence of total length on trophic position for walleye within the Inland Waterway .......................................................... 100 REFERENCES ........................................................................................................................... 102 vii LIST OF TABLES Table 1.1: Model set representing the multiple hypotheses evaluated to represent post-spawning walleye movement and demographics in the Inland Waterway, 2011-2013. Combinations of lake specific and time varying parameters for movement (φ), spawning-site fidelity (ψ), and fishing mortality (F) were evaluated using Deviance Information Criteria (DIC). F(.) is constant fishing mortality for each lake and time. The best fit model was also modified and fit without process error (*) to evaluate model support………………………………………………………………70 Table 1.2: Location specific post-spawning movement rates (with 95% credible intervals) estimated by models 1 and 3 under the base assumption for tag shedding rate (Ω = 0.14)……...71 Table 1.3: Fishing mortality rates (with 95% credible intervals) estimated from model 1 with and without process error using the base assumption of an instantaneous tag shedding rate of 0.14. Sensitivity of fishing mortality rates estimated using model 1 (best fit model) with three different instantaneous tag shedding rates (Ω = 0.04, 0.14, and 0.24)…………………………………….72 Table 1.4: Location specific post-spawning movement rates estimated from model 1 with and without process error using different assumed tag shedding rates (Ω = 0.04, 0.14, and 0.24)…..73 Table 1.5: Number marked and recovered by location and year for the 15 tag cohorts, which was used to inform the tag-recovery model. Cohorts 1-3 for each location correspond to the individuals tagged during spring-spawning in 2011-2013, respectively. Recovery years correspond to the annual fishing season that the fish were captured and returned. For example, cohort 1 from Burt Lake was tagged during spring-spawning in 2011 and the recovery years 1-3 correspond with the number of tagged individuals recovered and reported during the 2011-2013 fishing seasons…………………………………………………………………………………...79 Table 2.1: Location-specific values for parameters used to simulate movement and harvesting dynamics for Michigan’s Inland Waterway. The Abundance (SD) was calculated from a spring mark-recapture study in the Inland Waterway. The estimated fishing mortality, spawning-site fidelity, and post-spawn movement rates and their associated 95% credible intervals were taken from Herbst et al. (Chapter 1)……………………………………………………………………81 Table 2.2: Set of simulation scenarios used to evaluate the implications that movement patterns and differing exploitation rates have on walleye spawning populations within the Inland Waterway. Abundance, movement rates, and recruitment were location-specific and held consistent among all scenarios. Initial abundance was estimated for each location during a spring mark-recapture study in 2011. Natural mortality (M) was held constant at 0.3 and movement rates estimated in Chapter 1 were incorporated for all scenarios. The base scenario (1 and 8) represented perfect implementation of the current harvest control rule (i.e., constant total exploitation rate of 35%) with and without process error, respectively. In scenarios (2-7) movement rates and angling mortality rates (u-angling) were location-specific and taken from estimated posterior distributions (Chapter 1)…………………………………………………….82 viii Table 2.3: Example of the influence that post-spawn movement dynamics have on the total allowable catch (TAC) for each of the locations within the Inland Waterway. Spawning N refers to the estimated abundance during the spring-spawning assessment and the summer N is the abundance at each location after accounting for spearing harvest and post-spawn movements. The spearing TAC was determined using a spearing exploitation rate of 10% and the spawning abundance. The angling TAC was determined using an angling exploitation rate of 25% and the estimate of summer abundance. The angling TAC adjustments were determined by evaluating the difference between the two movement scenarios, with negative values indicating that fewer fish would need to be harvested to account for movement dynamics while maintaining the u = 0.35 for the spawning populations throughout the waterway……………………………………83 Table 2.4: Mean exploitation rates for each spawning population under different levels of uncertainty and spearing exploitation (us= 0.05, 0.10, and 0.20) and observed lake-specific angling mortality rates (ua). The uncertainty evaluated was process error on total instantaneous mortality (Z) over summer within the population dynamics model. (No error = process error excluded; Error = process error included; Diff = difference between the two approaches)……. 84 Table 3.1: Characteristics of lakes in the Inland Waterway. The percent littoral is defined as lake area ≤ 4.6 m total depth. All values are based off the 2013 mean. Total Phosphorus (TP) represents the bottom-to-surface mean. Crooked Lake Secchi and Chl-a are from 2011, which were the most recent data. Water chemistry data were provided by Tip of the Mitt Watershed Council (Kevin Cronk personal communication). Lake size and bathymetric information are unpublished MDNR data……………………………………………………………………………………………….89 Table 3.2: Catch-per-unit-effort (CPUE: fish/net) of fish in forage gill nets fished during MayAugust, 2011 and June-August, 2012 in the Inland Waterway. Only species that represent potential prey items for walleye and other predatory fish are listed, so incidental catches of predator species were omitted……………………………………………………………………90 Table 3.3: Seasonal diet composition (% wet weight (g)) for walleye from the Inland Waterway, 2011-2013. Seasons were defined as spring (April-June), summer (July-September), and fall/winter (October-March). N is the total number of stomachs that were analyzed for prey items. Unk. fish category consisted of unidentifiable bones and the “other” category consisted of miscellaneous fish species such as emerald and spottail shiners………………………………...91 Table 3.4: Stable isotope sample sizes (N), estimates of trophic position and percent littoral (SD) along with the standard ellipse area corrected for sample size (SEAc) for walleye within the Inland Waterway, in northern Michigan…………………………………………………………92 Table 3.5: Stable Isotope summary statistics (δ13C, δ15N, trophic position (TP), and percent littoral) for each species analyzed by lake. N indicates the total sample size collected, average length was measured in TL (nearest mm), values in parentheses represent the SD, covar is the covariance between carbon and nitrogen values, and * indicates composite samples…………..96 ix Table 3.6: Trophic position by length regression line specifics along with the overall mean length adjustment to determine biological influence associated with the significant interaction term that indicated trophic position was influenced by total length differently by location……………...101 x LIST OF FIGURES Figure 1.1: Map of northern Michigan’s Inland Waterway that consists of four lakes (Burt, Crooked, Mullett, and Pickerel) and four major connecting rivers (north to south through the lakes: Cheboygan River, Black River, Indian River, and Crooked River)………………………74 Figure 1.2: Conceptual model depicting the process how a single cohort (e.g., Burt Lake cohort 1) is tracked through time using the presented tag-recovery model. For example, after the initial tagging, which coincides with the spawning period, each individual within Burt cohort 1 has the ability to move to any location within the waterway or can remain in Burt Lake. Following that post-spawn movement the individuals then experience the population and observation processes that are representative of the location they moved to after spawning. Prior to time step t+1, individuals either exhibit spawning-site fidelity and return to their original tagging location (i.e., Burt Lake) or remain in the location they emigrated to. Following the spawning period those individuals once again have the ability to move freely throughout the waterway……………….75 Figure 1.3: Comparison of observed tag recoveries (red line) and the predicted distribution of tag recoveries with Bayesian p-values for each cohort during 2011-2013…………………………..76 Figure 1.4: Posterior (and prior) distributions of location specific fishing mortality rates (F) obtained from the best fit model with the base assumption for instantaneous tag shedding rate (Ω = 0.14)……………………………………………………………………………………………77 Figure 2.1: Map of northern Michigan’s Inland Waterway that consists of four lakes (Burt, Crooked, Mullett, and Pickerel) and four major connecting rivers (north to south through the lakes: Cheboygan River, Black River, Indian River, and Crooked River)………………………85 Figure 2.2: Conceptual diagram depicting the modeling process of how a spawning population (e.g., Burt Lake population) is tracked and projected through time in the Inland Waterway using the stochastic simulation model. Abundance is initially estimated using mark-recapture methods during spawning. Following the assessment of abundance ( ̂ ) the population is subjected to tribal harvest (i.e., spearing (us)) within the spawning grounds. After spawning, lake-specific spawning populations exhibit post-spawn movement (ϕ) and are subjected to total instantaneous mortality (Z = angling (ua) + natural (M)) in summer feeding locations. The fraction of the spawning population that survives during time t then returns to spawning grounds (i.e., spawning-site fidelity (ψ)) or remains in the location that they resided in during summer feeding. New additions represent immigrants from other spawning populations that do not exhibit ψ. During time t+1 the spawning populations are projected forward with the addition of constant recruitment that is specified at a level consistent with spawning population-specific harvest from time t. Locations abbreviations: BL = Burt Lake, ML = Mullett Lake, CL = Crooked Lake, PL = Pickerel Lake, and BR = Black River……………………………………………………………………………………………...86 Figure 2.3: Exploitation rates for each spawning population under the base scenarios where the harvest control rule was implemented without error (panel A; u = 0.35). Panel B represents the xi base scenario, but also includes a process error term on total instantaneous mortality (Z) within the population dynamics model. The red line indicates the target exploitation rate set forth as the waterway’s harvest control rule (u = 0.35)………………………………………………………87 Figure 2.4: Influence of varying levels of spearing exploitation rates (Scenario 2: us= 0.05 (low); Scenario 3: us = 0.10 (medium); Scenario 4: us = 0.20 (high)) and location-specific observed fishing mortality on exploitation rate (u) for each spawning population within the waterway. The exploitation rates under the base scenario that represents perfect implementation of the harvest control rule (HCR) is also provided for reference. The red line indicates the 95% CI for the maximum exploitation rate (u = 0.35), given the spawning population abundance estimate and standard deviation………………………………………………………………………………..88 Figure 3.1: Map of northern Michigan’s Inland Waterway that consists of four lakes (Burt, Crooked, Mullett, and Pickerel) and four major connecting rivers (north to south through the lakes: Cheboygan River, Black River, Indian River, and Crooked River)………………………93 Figure 3.2: Stable isotope bi-plots illustrating mean trophic position and percent littoral (± 1 SE) for multiple species within the communities of Mullett, Burt, Crooked, and Pickerel lakes within the Inland Waterway. Species-specific abbreviations include: ZOP = zooplankton, ALE = alewife, YEP = yellow perch, WAE = walleye, RGB = round goby, SPO = spottail shiner, CRAY = crayfish species, GST = gastropods, HEX= Hexagenia spp., and SMT = rainbow smelt……………………………………………………………………………………………...94 Figure 3.3: Stable isotope bi-plot representing the percent littoral and trophic position for walleye collected from each of the lakes within the Inland Waterway, MI. The lines enclose the standard ellipse areas (SEAc) and represent the total isotopic niche area for walleye within each of the lakes. The symbols are data points that indicate the isotopic signature of individual walleye from Burt (open triangle), Crooked (cross hatch), Mullett (open circle), and Pickerel (“x”) lakes……………………………………………………………………………………………...95 Figure 3.4: Plot of analysis of covariance results for the influence of total length on trophic position for each lake…………………………………………………………………………...100 xii CHAPTER 1 A State-Space Model for Estimating Walleye (Sander vitreus) Movement and Fishing Mortality in Interconnected Lakes Introduction Many aquatic ecosystems are open systems, allowing organisms to move among habitat types and geographic regions. Habitats are typically heterogeneous among subsystems, leading to demographic rates that vary across space. Given the prevalence of open systems and mobility of organisms, populations are often spatially structured and can exhibit diverse movement patterns. As a taxon, fish are among the more mobile inhabitants of aquatic systems, often demonstrating variable movement patterns and complex spatial structures that can complicate decisions related to harvest management and species conservation. Given these challenges, estimating movement rates among systems, and understanding the spatial structure of fish stocks has been an area of interest for ecologists and resource managers for decades (Hilborn 1990; Schwarz et al. 1993; Brownie et al. 1993; Schick et al. 2008; Hendrix et al. 2012; Molton et al. 2012, 2013; Li et al. 2014). Movement dynamics of fishes are frequently evaluated using mark-recapture and/or tagrecovery studies in which individuals are uniquely marked, released, and then later recaptured live or recovered via harvest (Hilborn 1990; Brownie et al. 1993; Schwarz et al. 1993; Pine et al. 2003). Multiple models are available to estimate movement and demographic rates from tagging studies. Commonly used approaches assume probabilistic movement, demographic, and recapture processes (e.g., Brownie et al. 1993; Schwarz et al. 1993), or deterministic movement and demographic processes with all stochasticity arising through the sampling process (e.g., 1 Hilborn 1990). One of the more widely used models was developed by Hilborn (1990), where the tag-recovery model embeds a biologically meaningful deterministic population model within a statistical estimation framework that uses a Poisson sampling model. More recently, extensions of the Hilborn tag-recovery model have been developed incorporating size selectivity (Anganuzzi et al. 1994), natural and fishing mortality (M and F) and tag shedding (Ω) (Aires-daSilva et al. 2009). Many extensions of the Hilborn model account for over-dispersion of the observed counts (i.e., var > mean) by using observation likelihoods such as the negative binomial or log-Poisson, thus accommodating greater uncertainty in tag-recovery data than allowed by the Poisson distribution (Dupuis 1995; Whitlock and McAllister 2009; Dupuis and Schwarz 2007; Hendrix et al. 2012). As a consequence, many applications of tag-recovery models in fisheries contain parameters relevant to the biology and management of fishes (e.g., M, F, Ω) but assume all variation in the data arises as a result of the sampling process, and variation arising from stochasticity in demographic processes is negligible. Often deterministic process models may not be plausible approximations, given that vital rates for both individual animals and populations can exhibit considerable variation through space and time (Ogle 2009; Hansen et al. 2011; Bjorkvoll et al. 2012). Therefore, it is important to incorporate stochasticity in the underlying population model, and inclusion of both process and observation uncertainty will increase the realism of tag-recovery model applications in fisheries. A state-space model is a special class of a hierarchical statistical model for time series data that provides a rigorous approach for modeling stochastic biological and observation processes (Schnute 1994; King 2014). This framework also provides a flexible approach for tailoring biological process models to the life history of a study organism (Thomas et al. 2005). State-space models have been used to estimate demographic and movement parameters in mark- 2 recapture studies (e.g., Gimenez et al. 2007; Kéry and Schaub 2011), but have seen less application for estimating movement parameters of spatially-structured fish populations using tag-recovery data (e.g., extensions of the Hilborn model). Applications of state-space models in fisheries commonly employ Bayesian estimation and model fitting techniques, providing additional flexibility by allowing one to easily constrain parameter values over realistic ranges or explicitly incorporate prior data recorded from other time periods, populations, or species (Whitlock and McAllister 2009). Use of prior information in conjunction with Bayesian estimation techniques has gained popularity in ecological studies because estimating parameters from complex state-space models still remains technically challenging (Schnute 1994; Kéry and Schaub 2011; King 2014). Moreover, this approach is advantageous because it allows sharing of information across data sets when specific pieces of information are lacking on the population or site of interest (Whitlock and McAllister 2009; Kéry and Schaub 2011; Hendrix et al. 2012), a common scenario for agencies tasked with managing populations over vast geographic areas with limited resources. Despite the strengths of the state-space model framework, many applications have not yet incorporated important aspects of fish behavioral ecology. Many fish species exhibit regular seasonal or inter-annual movement patterns associated with reproductive events and movement to seasonal feeding habitats, especially in open systems. Movements associated with spawning events are often nonrandom and spawning-site fidelity is a common life-history attribute for a wide variety of fish species (Moyle and Cech 2004). For example, the walleye (Sander vitreus) is a mobile species that typically exhibits seasonal movements from spawning to feeding areas. However, these movement patterns vary among systems in the extent of directed movement displayed (Rasmussen et al. 2002; DePhilip et al. 2005; Weeks and Hansen 2009). Although walleye post-spawn movement appears to be context 3 dependent, individuals are regularly captured in the same general location during the annual spawning period, which suggests that walleye exhibit some degree of spawning-site fidelity (Crowe 1962; Olson and Scidmore 1962). In general, the structure of current tag-recovery models does not incorporate explicit across-year returns to a specific location or within year movement among locations, and thus ignores these life history characteristics (e.g., spawningsite fidelity) that may be important for understanding movement dynamics and spatial structure of walleye populations. Despite the importance of walleye as a game species, relatively few studies have quantified their movement rates (Rasmussen et al. 2002; Weeks and Hansen 2009; Vandergoot and Brenden 2014), which is likely due to logistical challenges as well as the limitations in analytical tools to account for these complicated movement patterns. The goal of this study was to understand and quantify the movement dynamics of walleye in a set of large interconnected lakes and river systems in northern Michigan. Specific objectives of this study were to: 1) develop a tag-recovery model that accounts for the biology of our study system and integrates prior sources of data to estimate movement and demographic parameters and 2) quantify the movement rates of a walleye chain-lake system in northern Michigan during 2011-2013. To accomplish these objectives we developed a state-space tag-recovery model that adapts the general framework of Hilborn (1990), described further by Quinn and Deriso (1999), to account for important biological characteristics of this system (e.g., movement dynamics, spawning-site fidelity), while integrating prior data sources that allowed us to estimate important demographic and fishery parameters (e.g., exploitation rate) in each lake. This model was implemented in a Bayesian estimation and inferential framework, providing a flexible approach for understanding 4 dynamics and permitted stochasticity in both underlying biological processes and in observation processes generating the tag-recovery data (Gimenez et al. 2007). Methods Study area Michigan’s Inland Waterway is an interconnected chain of lakes located in the northern Lower Peninsula that consists of four large lakes (Burt, Crooked, Mullett, and Pickerel) interconnected by a series of rivers and smaller tributaries (Figure 1.1). The Cheboygan Lock and Dam on the Cheboygan River and the Alverno Dam on the Black River located at the northern portion of the Inland Waterway restrict fish passage and are considered closed to emigration towards Lake Huron or further upstream within the Black River (Figure 1.1). The lakes and rivers of the waterway are oligotrophic, provide various levels of suitable walleye spawning substrate and prey resources, and range from 4.4 km2 (Pickerel Lake) to 70.4 km2 (Burt Lake) in total size (Hanchin et al. 2005a; Hanchin et al. 2005b). The Inland Waterway was separated into five spatial strata consisting of the four lakes and the Black River, for the purpose of this study. Boundaries of the spatial strata were defined as 1) the Black River, 2) Mullett Lake including the Cheboygan River, 3) Burt Lake including Burt Lake, Indian River, Sturgeon River, and the Crooked River, 4) Crooked Lake including Crooked Lake and the Crooked-Pickerel narrows to the mid-point between Crooked and Pickerel lakes, and 5) Pickerel Lake including Pickerel Lake and the other half of the Crooked-Pickerel narrows nearest to Pickerel Lake. The divisions of these waterbodies into the specific strata were based on the four lakes, and the rivers were categorized based on proximity to a specific lake and biological information gained from past walleye studies in the Inland Waterway and input from 5 local biologists (Michigan Department of Natural Resources, unpublished data). For example, the Cheboygan River was categorized into the Mullett Lake strata because the majority of walleye captured in the river during spring sampling were collected within 150m of Mullett Lake. Tagging and recovery data Adult walleye, defined by expression of gametes or total length ≥ 381mm, were captured in the spring (mid-March to early-May) during the walleye spawning season using electrofishing, fyke nets, and trap nets throughout the Inland Waterway, 2011-2013. Following capture, walleye were marked with individually numbered, size 12 jaw tags that were affixed to their upper mandible. Tags also were labeled with a mailing address for return, and approximately half of the jaw tags affixed were $10US reward tags to increase reporting rate. Information recorded for each individual during tagging included location, date of initial marking, and total length (mm), and sex if gametes could be expressed. Tag recovery data were provided to the Michigan Department of Natural Resources through a voluntary angler tag return program during the 2011-12, 2012-13, and 2013-14 angling seasons. The information collected from each tag recovery included date and location of capture. In addition to the monetary reward, project collaborators advertised the return program to the angling community through public outreach events, press releases, and signage at access points to encourage tag returns. 6 Model Structure General Approach Application of hierarchical statistical models in applied ecology over the last decade has resulted in substantial advancements in the complexity and realism of ecological data analyses (Ellison 2004; Clark 2007; Royle and Dorazio 2008). Hierarchical approaches facilitate the integration of stochastic process models from ecological theory with real datasets via explicit modeling of both observation and ecological processes (Royle and Dorazio 2008; Kéry and Schaub 2011). Royle and Dorazio (2008) described hierarchical models as the joint probability distribution of conditional probability models representing observational and ecological processes: ( ( ) ) , for observed data y, partially observed latent state variable X (the true quantity of interest), and potentially vector-valued parameters governing the observation ( ( ) and ecological processes ). Bayesian statistical tools are commonly used to combine these probability models to make inferences on model parameters via the joint posterior distribution: ( where ( ) ) ( ( ) ( ), ) is the joint prior distribution of all observation and process parameters (i.e., parameters and hyper-parameters). In the context of modeling fish movement among spatial strata, ( ) represents the stochastic process that determines how many individuals are available to be caught during an angling season in a given geographic strata, which is governed by the movement and demographic parameters. In contrast, ( ) represents the probability distribution for observing y tag recoveries from a given tagging cohort given the true 7 number of fish available for harvest (X) and angling harvest and tag-reporting processes. We developed a state-space tag-recovery model, in which the population states (X) are linked through time in a Markovian fashion (e.g., Schnute 1994; King 2014). Thus, the spatial distribution of fish in each lake available for harvest at time t+1 depends on the distribution and abundance of fish at time t and the movement and demographic parameters, and specific realizations of process errors (Figure 1.2). The tag-recovery state-space model parameters that are associated with the process and observation models and their descriptions are provided in Table 1.1. Population process model The process component of our state-space model governed the spatial-temporal dynamics of movement and survival of fish from each tagging cohort. Specifically, the number of fish from each unique release group (i.e., cohort) available for harvest on summer feeding grounds in a given year was a latent variable (X). Changes in X were modeled as a function of the number and spatial distribution of fish from that group at the previous time step and the parameters driving demographic processes of movement and apparent survival. These processes were governed by the following general model, in which fishing mortality and movement rates are year specific: ( ∑ ) (∑ ( ) ) . Here Xj,l,i,t represents the number of fish from tagging cohort j released on spawning grounds at site l that are present and available for harvest on summer feeding grounds at site i during year t. 8 In this study there were three release cohorts (j = 1,..,3) at each of 5 spatial strata (i = 1,…,5) resulting in 15 unique release groups, and three harvest recovery years (t = 1,…,3). Moreover, Rj,l,t represents the number of tagged fish released in cohort j at spawning site l at the start of year t, and is the apparent annual survival rate for walleye at site l during time t. Also note that the number of released fish ( (i.e., when ), whereas this term is only non-zero when equals the release year for a given tag cohort (i.e., when ) equals zero for tag cohorts not released during recovery year t ). Thus, the term is only necessary in the above state equation because all fish are not tagged and released during the first recovery year of the study. We also evaluated simpler models representing alternative hypotheses where parameter values were constrained to be equal across space and/or time. Our process model assumes that all mortality occurs after fish move to summer feeding areas. The process model also assumes that all fish in a given feeding area during recovery year experience the same conditions and thus are experiencing the same apparent survival. Similarly, fishing mortality is operating at the site level during the summer, where fish in the same site are exposed to similar levels of fishing mortality, regardless of their unique release group. Because the processes governing movement and survival dynamics (e.g., Hilborn et al. 1990; Hendrix et al. 2012) are unlikely to be deterministic, we incorporated a multiplicative process error that represents the cumulative result of stochastic variation in all mortality and tag-loss processes. Specifically, we assumed process error was acting on total instantaneous mortality in a manner that was lake and time specific: ( where M = 0.3, Ω = 0.1375, ( ) ), and , . Natural mortality (M) was assumed and held constant at a value consistent with estimates of M from walleye populations in northern 9 Wisconsin (Hansen et al. 2011). Our base assumption for the tag shedding rate (Ω) was reflective of an estimate from a walleye mark-recapture study conducted within our study area in 2001 (Hanchin et al. 2005a, 2005b). This model formulation treats tag-loss as a component of instantaneous total mortality of the tagged population, and as such Z does not represent true mortality but apparent total mortality. Thus our model has no way of separating out components of process error related to tag-loss and true mortality, as our data do not permit such partitioning. Examining the above state equations governing the distribution and abundance of tagged fish from each release group aids the interpretation of model dynamics. Thus, ∑ represents fish from tag cohort that survived at their initial release location during time t-1, plus the sum of all fish that survived at other sites during time t-1 and then returned to spawn at their initial release location l at time t. This sum therefore represents the number of fish that will be available to move from their initial spawning locations (i.e., fish exhibiting spawning-site fidelity) at time t to summer feeding grounds at site i, where is the proportion of the fish that will make this movement. Moreover, ∑ ( ) represents the sum of fish from tag cohort j that survived at sites other than their release location during time step t-1, and subsequently remained at these locations for spawning at time t (i.e., failed to exhibit spawning-site fidelity). This sum therefore represents the number of fish that will be available to move from their current spawning location (not their initial spawning location) to summer feeding grounds at site i, where is the proportion of the fish that will make this movement. Thus fish that do not exhibit spawning-site fidelity remain to spawn at the 10 location where they summered during the previous time step, and the overall model represents the number of fish from each release group that are present and available for harvest on summer grounds at site i during recovery year t. Also note that all of our models assumed spawning-site fidelity was constant through time, whereas movement rates from spawning grounds to summering grounds were allowed to vary through time for some models. Tag-recovery observation model While the stochastic process model above drives movement and survival dynamics of fish from each tagging cohort, the observation model describes the processes generating observed tagrecovery data. We assume tag-recovery is a stochastic process, where the number of tags recovered from each release cohort at each site and time is conditional on the number of fish present with tags and the parameters driving harvest and tag reporting: ( where ) represents the number of walleye tags recovered at site i during time t from fish released in tag group j at site l. The mean of the Poisson distribution for tag recoveries was determined by the number of fish available for harvest, the annual exploitation rate, and the tagreporting rate: , where ( ) is the annual exploitation rate for walleye at site i during time t. Because this model is assuming recoveries are coming from summer feeding grounds we assume here that all fish present in a given space-time combination are experiencing the same exploitation rate, regardless of which tag cohort they belong to or where they spawn. Reporting 11 rate ( ) was assumed constant over space and time and was estimated using auxiliary reward tag data (see below). Prior distributions for model parameters We used existing data to develop informative prior distributions for model parameters where available, and used a diffuse prior distribution for the parameters that were of primary interest for this analysis. We used pooled catch-at-age data from walleye collected throughout the Inland Waterway in 2011 to develop an informative prior for fishing mortality using results from a catch-curve analysis. We loge transformed the catch curve equation to estimate instantaneous total mortality rate (Z) using linear regression (Quinn and Deriso 1999). From the catch-curve analysis ̂ was the maximum-likelihood estimate of instantaneous total mortality, which has an asymptotically normal sampling distribution ( ̂ ( ̂ ) ). We assumed that ) and thus ̂ natural mortality was constant over the catch-curve study period ( . Since linear functions of normal random variables are themselves normally distributed (Rice 2007) we used results from catch-curve analyses to derive an informative normal prior for linear function of the normally distributed random variable ̂ (Appendix 1.1): ( ). To avoid impossible or unrealistic draws from the prior for (MCMC) sampling discarded any samples of , Markov Chain Monte Carlo ≤ 0 and ≥ 5. We lacked prior information about the magnitude of process errors, so we assumed a uniform prior over a restricted range: 12 as a Although this prior is uniform, the bounds of the uniform distribution can be thought of as informative in this case because we are restricting the process error values over a relatively small numerical range. While this range is numerically restrictive, it contains all biologically plausible values of the process error; process error approaching produces u-shaped distributions of apparent annual survival ( ) where nearly all individuals in the tag group survive or die (or shed tags) each year. This is a biologically unrealistic situation for walleye in northern Michigan, thus the uniform prior used constrains the process error standard deviation to biologically reasonable values while at the same time reflecting ignorance over the values of within a plausible range. We used data from live recaptures during annual tagging operations to develop informative priors for spawning-site fidelity parameters ( ). Specifically, the number of fish tagged on spawning grounds that were recaptured at their initial release site in subsequent years, was treated as a Binomial random variable with success probability for site . The conjugate prior for a Binomial parameter is a Beta distribution, and the Uniform distribution represents a special case ( ( ) ). Moreover, using an uninformative prior for a Binomial parameter results in a closed-form Beta posterior distribution for the Binomial probability ( ( )), where x = number of successes from n Bernoulli trials. Thus, we used this approach to turn the proportion of tagged fish recaptured on their original spawning release area into an informative Beta prior ( )) for the spawning-site fidelity parameter for a given site ( ( ), where n represented the number of fish tagged from spawning site recaptured on any of the spawning grounds during tagging operations for subsequent spawning seasons, and x represented the number of these fish recaptured at their original spawning ground release locations. For example, 485 walleye released on spawning grounds in Burt Lake were recaptured during 13 tagging operations in subsequent spawning seasons, and 479 of these fish were recaptured within Burt Lake. This resulted in a ( ) prior for ( ( )). This approach was used to turn the posterior distributions from Bayesian estimation of site-fidelity parameters using live-recapture data into informative priors for spawning-site fidelity for all sites when fitting the full state-space model: ( ), ( ( ), ( ( ), ), ). We developed an informative prior distribution for reporting rate using data collected during high-reward walleye tagging studies conducted in Crooked, Pickerel, and Burt Lakes in 2001 and within the entire Inland Waterway in 2011 (MDNR unpublished data). The standard tag return rate (λ) and its variance were estimated via the ratio of the recovery rate of standard tags to the recovery rate of high-reward tags; these methods are described further within Henny and Burnham (1976), Conroy and Blandin (1984), and Pollock et al. (1991). The estimate (mu) and variance of λ were then used to develop an informative Beta prior for the reporting rate parameter: ( ) ( ( ) ). Because we lacked prior information on movement from spawning to feeding grounds among lakes and because these were our primary targets of inference, we used diffuse priors for all parameters. Two sets of constraints must be met for the vector of movement rates away from spawning site at time : 1) movement rates away from a site must be bound on the interval 14 [0,1], and 2) all movement rates leaving site at time must sum to one. For the vector of movement rates out of a given site at time t we used a diffuse Dirichlet distribution, which is a multivariate generalization of the Beta distribution that fulfills the necessary set of parameter constraints (Gelman et al. 2004). Thus, we specified a vague Dirichlet prior for each ( where ), = 1 for all sites receiving fish from site l at time t. This effectively allocates individuals uniformly across all receiving sites at time t (Royle and Dorazio 2008). To implement this prior we simulated independent Gamma(1,1) random variables, and expressed movement rates out of site l as functions of these random variables (Royle and Dorazio 2008): ( ) ⁄∑ . Model set We developed a set of 8 models representing hypotheses of how movement (φ) and fishing mortality (F) rates potentially vary by location and time. In particular, our model set allowed for both site and time specific movement and fishing mortality parameters, but all models assumed spawning-site fidelity was lake-specific and constant over time (Table 1.2). To evaluate relative support for our alternative models we used deviance information criteria (DIC; Spiegelhalter et al. 2002), which is calculated as a function of the posterior distribution of model deviance and the number of effective parameters (pD). 15 Model fitting and evaluation Models were fit using OpenBUGS (Bayesian inference using Gibbs sampling) software (http://www.openbugs.net) called from the R2OpenBUGS package within R (R Development Core Team 2010). Samples from the posterior distributions of all model parameters were generated using Gibbs sampling, and all analyses used three Markov Chain Monte Carlo (MCMC) chains with random starting values for model parameters. Preliminary analyses suggested all MCMC samplers converged to the posterior distributions after approximately 100,000 iterations. Thus, for each chain we used a burn-in period of 150,000 iterations that were discarded followed by 200,000 samples that were retained, resulting in posterior distributions described by 600,000 samples for each model parameter. All chains were evaluated for convergence and mixing using the Gelman-Rubin statistic (Gelman and Rubin 1992) and visual inspection of traceplots and posterior density plots for all model parameters. We evaluated fit of the top state-space model to our tag-recovery data using Bayesian p values, which provided comparison of the posterior predictive distributions of model predicted quantities with the observed tag-recovery data (Meng 1994). Specifically, we calculated a Bayesian p value for the omnibus chi-square statistic (Gelman et al. 2004), where the posterior predictive distribution of the chi-square statistic was a weighted measure of discrepancy between the predicted and observed number of total tag returns from all sites and cohorts over all posterior samples of model parameters. Bayesian p values close to 0.5 represent a good fit of the model to the data, since on average the predicted values are less than or greater than the observed value with equal frequency (Whitlock and McAllister 2009). While the omnibus chi-square statistic is a measure of fit over the entire model, we were also interested in evaluating fit of our model to tag-return data from each tagging cohort. Thus, we calculated the posterior predictive 16 distributions for the sum of all tag returns across all sites from each specific release group and compared this to the observed tag returns using Bayesian p values. This provided an indication of specific areas where model assumptions may have been violated or areas where the model simply did not predict the raw data well. Results Model Selection Eight different models were fit using the walleye tag-recovery data to evaluate support for hypotheses that represented various combinations of how movement (φ) and fishing mortality (F) varied by location and time. The top model as indicated by DIC values included spawningsite fidelity, movement, and fishing mortality rates that were location specific and constant during the three year study (Table 1.2). The second best model (i.e., model 3) showed limited support (ΔDIC = 6.5) for location specific spawning-site fidelity and movement rates, and fishing mortality rates that varied by location and year. We evaluated model fit and the posterior distributions of parameters obtained from the top model because of its DIC value. However, estimated movement and fishing mortality rates obtained from the top two models resulted in very similar parameter estimates (Table 1.3). Evaluation of the top model showed good model fit to observed walleye tag returns (χ2 = 0.83, P = 0.56). Furthermore, fit of the model to tag-return data for each of the 15 cohorts demonstrated that the posterior predicted distributions for tag recoveries fit well with the actual numbers observed for nearly all tagging cohorts (Figure 1.3). The few exceptions were cohorts that had lower numbers of observed tag recoveries (i.e., Mullett Lake cohorts 1 and 3, and Black River cohort 3), which had p-values that strayed away from the optimal value of 0.5 (Figure 1.3 and Appendix 1.2). 17 Demographic parameters Movement rates Walleye within the Inland Waterway exhibited asymmetrical post-spawning movement patterns, based on estimates obtained from the best fit model. Moreover, estimated post-spawn walleye movement rates were consistent across both of the top two models (Table 1.3), thus we only describe results from the top model. The Black River had the highest post-spawning emigration rate. Of the cohorts initially tagged in the Black River approximately 82% (95% CrI: 0.52 - 0.96) emigrated to other areas for summer feeding (Table 1.3). Of the 82% exiting the Black River after spawning, the majority (Mean = 77%; 95% CrI: 0.48 - 0.92) moved into Mullett Lake (Table 1.3). Among the lake sites, post-spawning emigration was highest from Pickerel Lake, where an estimated 30% (95% CrI: 0.61 - 0.80) of the population left the lake. Bi-directional post-spawn movement of walleye between Crooked and Pickerel lakes occurred more frequently than other combinations of locations. Post-spawn movement of walleye from Crooked Lake to Pickerel Lake was less substantial (Mean = 6%; 95% CrI: 0.03 - 0.10), but 19% (95% CrI: 0.13 0.25) of fish spawning in Pickerel Lake moved to Crooked Lake during the feeding season (Table 1.3).Walleye cohorts initially tagged in Burt and Mullett lakes had greater overall site fidelity, with 90% (95% CrI: 0.87 - 0.93) and 86% (95% CrI: 0.63 - 0.94) remaining in those locations throughout the year, respectively (Table 1.3). Fishing mortality The number of fish tagged and number of returns varied widely between locations in the watershed, and as such, the level of information provided for parameter estimation varied. A comparison of the difference between the prior and posterior distribution for fishing mortality 18 (F) for each location (Figure 1.4) indicated that the tag recovery data were informative for estimating location-specific Fs for most sites. Estimated fishing mortality rates fell into two broad groups; the Black River, Burt Lake and Mullett Lake all had an estimated F between 0.15 and 0.19 (Table 1.4), whereas F estimates in Pickerel and Crooked Lakes were 0.28 and 0.31(Table 1.4), respectively. The posterior distributions of F for most lakes were symmetrical and reasonably narrow (Figure 1.4). However, the posterior distribution of F in the Black River was asymmetrical and multimodal (Figure 1.4), with a 95% credible interval ranging from 0.02 to 0.28 (Table 1.4), suggesting that the low number of tag returns from the Black River resulted in only partial identifiability for fishing mortality at that site. The mean year-specific F values by location provided by model 3 were generally similar to F values estimated using model 1, which provided a single F for the study period for each location. The Black River was again the exception, where the estimated fishing mortality rate using model 1 was 0.15 and the mean of the annual estimates from 2011-2013 using model 3 was higher (F = 0.22). Issues with estimability and inconsistencies of fishing mortality estimates from Black River were likely due to the low number of tag returns from this location (Appendix 1.2). For other locations, there was limited annual variability in estimates of F obtained from model 3, providing further support for the best fit model, which assumed constant F over time. A single exception occurred in Pickerel Lake in which the fishing mortality rate estimates differed by 0.10 from 2012 (F = 0.30) to 2013 (F= 0.20). Sensitivity to process error and tag shedding rate assumptions The tag-recovery data used to inform the model allowed for the estimation of the process error that was incorporated into the population dynamics and observation models. The posterior distribution for the process error parameter was symmetric with a mean of 0.76 (95% CrI: 0.44 – 19 1.43), which differed from the uniform (0,3) distribution that was used as the prior distribution. In addition to directly representing variability in demographic processes, the process error term was also responsive to other parameters such as tag shedding rates. For example, the process error increased from a mean of 0.76 to 0.94 with the decrease in the tag shedding rate from 14% to 4%, while estimates of movement and fishing mortality rates remained consistent during the same modification (Table 1.4 and 1.5). Thus, inclusion of process error appeared to facilitate more robust estimation of movement and fishing mortality parameters than the model without process error. Moreover, overall support substantially declined after modifying the structure of the best fit model to exclude process error (ΔDIC = 112; Table 1.2), indicating the added value of including process stochasticity in the model structure. The best fit model was robust to varying estimates of instantaneous tag shedding rates. Fishing mortality rate estimates varied by less than 0.01 in response to increasing tag shedding rates from 0.04 to 0.24 (Table 1.4). Variation in estimated movement rates was generally low (Table 1.5) in response to this range of tag shedding rates. The process error parameter was influenced more by the variation in instantaneous tag shedding rates, increasing when the value for tag shedding rates (Ω) decreased. The estimated process error standard deviation when using the low Ω value was 0.94 (95% CrI: 0.53 - 1.88) and at the high Ω was 0.65 (95% CrI: 0.37 1.17), illustrating the variation in process error following a change in tag shedding rate from 4% to 24%. In addition, when process error was removed from the model structure the parameter estimates of F were more sensitive to variable levels of Ω (Table 1.4). Movement rates did not exhibit a similar pattern, remaining the same after the exclusion of process error from the population-process model structure. The only exception was the Black River, which varied slightly (< 2%) and was likely associated with a lower number of tag-returns (Appendix 1.2). 20 Discussion This study expanded upon previous extensions of the commonly used Hilborn (1990) tagrecovery model by developing a state-space formulation to accommodate spawning-site fidelity, a common life history trait of many fish species. The state-space tag-recovery model structure developed accommodates temporal and spatial variation in demographic and movement rates (i.e., F and φ), and includes process stochasticity to help alleviate inferential sensitivity associated with commonly used but incorrect assumptions like constant and known rates of natural mortality and tag-shedding. In addition, the state-space framework presented here benefited from the use of Bayesian estimation techniques, which provided the model the flexibility to incorporate site- and species-specific knowledge through the use of prior distributions while estimating population rates of interest such as movement (φ) and fishing mortality (F). The Bayesian approach also facilitated inclusion of prior information while accounting for uncertainty in that knowledge (i.e., through prior distributions) and thus we avoided simply assuming fixed parameter values for quantities not likely to be estimable using only the tag-recovery data (e.g., spawning-site fidelity). Thus we were able to embed more realistic biological dynamics into the model structure while using existing auxiliary information to aid model fitting, which are recognized benefits of Bayesian implementations of state-space models (Buckland et al. 2000; Buckland et al. 2007). Furthermore, this approach was complemented by formal statistical evaluation of hypotheses of the spatial and temporal structure of model parameters (e.g., constant versus time-varying φ and F) using accessible Bayesian model selection approaches, thus making the general approach useful under a wide range of biologically plausible conditions within both freshwater and marine environments. 21 The state-space tag-recovery model presented in this study allowed stochasticity to be incorporated into both the population dynamics model and the observational sampling process. This is different from other commonly used extensions of the Hilborn (1990) model that assume all uncertainty is within the observation processes of tag-recovery, where additional uncertainty is typically accounted for by assuming an over-dispersed likelihood (Anganuzzi et al. 1994; Aires-da-Silva et al. 2009; Hendrix et al. 2012). For example, Hendrix et al. (2012) recommended that a negative binomial or a log-Poisson likelihood approach be used to account for the statistical uncertainties associated with individual-level heterogeneity in recapture probabilities and group movement patterns. Our approach, however, accounts for additional uncertainties by incorporating stochastic variation in both the population dynamics and observational processes, providing a more realistic framework for developing robust estimates of movement parameters, additionally, our approach is flexible enough to be adapted to specific behaviors (e.g., spawning-site fidelity). Moreover, inclusion of process error allows for stochastic uncertainty in realizations of the underlying numbers of fish available for harvest at each space-time stratum (which in turn drive expected tag returns), and thus avoids assumptions of strictly deterministic (and correctly specified) dynamics when in fact demographic processes are often highly variable (Schmalz et al. 2011). Thus, including process error into the structure of tag-recovery models allows for variability in the behavior of individual fish or spatialtemporal variation in the demographic rates that are not explicitly contained within the population dynamics model. Accounting for both sources of uncertainty is thus advantageous over commonly used historical approaches that relied on deterministic population models with fixed dynamics (Buckland et al. 2000). Incorporating uncertainty into both the process and observation models of the state-space tag-recovery model allows for the stochasticity in the 22 demographic processes to be absorbed into the process error, which can increase the robustness of parameter estimates (i.e., φ and F observed in this study). Estimates of demographic rates from fish populations can be biased because of uncertainty in the magnitude of tag shedding (Isermann and Knight 2005; Aires-da-Silva et al. 2009; Koenigs et al. 2013). For example, previous studies have generally shown that estimates of movement and fishing mortality rates are sensitive to tag shedding (Isermann and Knight 2005; Aires-da-Silva et al. 2009). Immediate or short-term tag shedding is often low for walleye (< 0.05%), but long-term tag shedding for walleye is more variable and has been estimated to range between approximately 5 and 50% annually (Hanchin et al. 2005; Isermann and Knight 2005; Koenigs et al. 2013; Vandergoot et al. 2012). Given the high uncertainty in tag shedding rates, we assessed a gradient of instantaneous tag shedding rates (Ω = 0.04-0.24) and demonstrated that our movement rate estimates were robust to these assumed values. In fact, the difference in postspawn movement rates was ≤ 2% for each of the tag shedding scenarios evaluated. The insensitivity of movement rates to variable levels of tag shedding was unexpected based on results from previously cited tag-recovery studies. For example, Aires-da-Silva (2008) reported that estimates of mean movement rates for blue sharks were highly sensitive to high levels (Ω = 0.22) of tag shedding when compared with estimates associated with lower tag shedding rates (Ω ≤ 0.11), with movement varying as much as 0.14 under the different assumptions of tag shedding rates. The demographic rates of interest were robust within our best fit model, which is likely the result of our model structure allowing for additional stochasticity in the instantaneous total mortality through the inclusion of process error instead of assuming total mortality is a function of tag-loss assumptions within a rigid deterministic model governing movement and demographic dynamics. 23 Walleye Dynamics in the Inland Waterway Walleye movement has been shown to vary widely among systems studied. For example, Rasmussen et al. (2002) found that at least half of all walleye present at spawning could emigrate to another site within one week in a chain-lake system, whereas Weeks and Hansen (2009) found that the majority (82%) of walleye tagged were recaptured in the same lake. Our results demonstrated that walleye within the Inland Waterway exhibited asymmetrical post-spawning movement patterns similar to walleye in other lake-chain systems (Rasmussen et al. 2002). Although our study was not designed to determine factors governing movement rates, we expect that walleye populations in lakes with suitable spawning substrate and abundant prey resources would not benefit from migrating great distances to spawn and/or feed. Alternatively, if spawning substrate and adequate forage are spatially separated it would be advantageous for those walleye to migrate greater distances in search of quality habitats, thereby increasing chances of juvenile survival and/or adult growth. Fishing mortality rates varied within the waterway, but were consistent over the three year study and were within the range reported for other walleye populations (Schmalz et al. 2011). Within the Inland Waterway, Crooked and Pickerel lakes had the highest estimated fishing mortality rates (F = 0.31 and 0.28 respectively); however, neither of these rates exceeded 35%, which is commonly viewed as an upper limit reference point for safe harvest of walleye (Schmalz et al. 2011). Estimated fishing mortality rates in the other lakes and in the Black River were in the range of 0.15 to 0.31, suggesting that exploitation is not the primary factor limiting abundance of adult walleye in these systems. We observed little inter-annual variability in estimated fishing mortality rates from the top two models identified in this study. This is in 24 contrast to estimates reported by Hansen et al. (2011) for walleye in Escanaba Lake, WI where the annual exploitation on age-3 and older walleye differed as much as 50% annually, demonstrating that annual variability for this mortality source can be high. In summary, this study expanded a commonly used tag-recovery modeling framework to incorporate spawning-site fidelity and additional uncertainty associated with the population dynamics processes into the model structure using a state-space framework. We used Bayesian estimation techniques to facilitate inclusion of existing information while accounting for uncertainty through the use of prior distributions. We determined that post-spawn walleye movement patterns in the Inland Waterway were spatially asymmetrical but were consistent over the study period. Furthermore, our movement and fishing mortality estimates were robust to changes in assumed rates of tag loss. Given the prevalence of open systems and organisms with complex life-history behaviors, flexible modeling frameworks that incorporate stochastic process dynamics and are readily adaptable to different species and systems are important additions to approaches commonly used to model tag-recovery data in fisheries. State-space models like the one presented in this study thus provide a state-of-the art framework that will permit scientists to robustly estimate demographic parameters governing movements and mortality for mobile species (King 2014), which will ultimately provide rigorously evaluated information to aid management decisions for spatially structured fish populations. 25 CHAPTER 2 The Influence of Movement Dynamics on Harvest Allocation for Walleye (Sander vitreus) in Intermixed Fisheries in a Chain of Lakes Introduction Fish are highly mobile organisms, often with complex movement patterns that result in spatially structured populations. Difficulties associated with delineating a population’s spatial structure and accounting for intermixing in assessments presents several challenges in the fishery management process. Two prominent issues associated with managing spatially structured populations include biased estimates of abundance or other demographic parameters and difficulties in implementing harvest control rules. It is imperative, however, to address these issues because ignoring spatial structure in fisheries assessments can lead to overexploitation (Ying et al. 2011; Molton et al. 2012, 2013; Guan et al. 2013) and has led to the collapse of commercially harvested species (Morishima and Henry 1999; Ames 2004; Fu and Fanning 2004; Hutchinson 2008). The risk of stock collapse illustrates the significance of understanding and addressing spatial structure when making management decisions. Lake chains within the ceded territory of the northern Great Lakes region support multiple fisheries and provide examples of complex systems where managers are faced with spatially structured fish populations (Rasmussen et al. 2002; Chapter 1). The Inland Waterway, for example, is a lake chain in northern Michigan that supports walleye populations that exhibit asymmetrical post-spawn movement patterns (Chapter 1). The management strategy within this waterway is based on an agreement (2007 Inland Consent Decree) between state and tribal resource managers that treats each lake within the waterway as a closed system, consisting of a discrete fish stock. For each lake, harvest is limited using a harvest control rule of a constant 26 total annual exploitation rate (u = 0.35) that is allocated to tribal exploitation (us = 0.10) and recreational angling exploitation (ua = 0.25). The premise of this strategy is that an exploitation rate of 0.35 poses a low risk of overexploitation, as summarized by Schmalz et al. (2011). The total allowable catch (TAC) for each location within the waterway is set using the harvest control rule and is based on lake-specific mark-recapture estimates of adult abundance, or if markrecapture estimates are not available, an estimate of abundance derived using regression-based estimates (Hansen 1989; Nate et al. 2000). Allocation and monitoring of harvest is challenging because of seasonal intermixing among lakes (Chapter 1) and the presence of two fisheries that occur over different temporal scales. Tribal spearing harvest occurs during the spring spawning season (late-March through April) when fish are in their natal locations, whereas recreational anglers harvest mixed populations during the state-regulated fishing season from May through mid-March. In addition, walleye populations within the Inland Waterway have heterogeneous levels of recruitment and rates of growth (Hanchin et al. 2005a; Hanchin et al. 2005b; Michigan Department of Natural Resources unpublished data). The combination of these factors creates challenges for managers when developing and implementing management strategies. Management decisions are typically based upon assessment methods that estimate system state values (e.g., abundance) from spatial areas that are assumed to contain a single spawning stock (Quinn and Deriso 1999). However, when movement causes intermixing and when individuals from multiple spawning stocks are present within an assessment area, typical singlestock assessments can lead to over estimation of abundance, which in turn can lead to overexploitation (Ying et al. 2011; Molton et al. 2012, 2013; Guan et al. 2013). Movement rates are commonly unknown or difficult to quantify, and as such, assessment models often ignore or make simple assumptions about movement rates (Booth 2000; Walters and Martell 2002; Cope 27 and Punt 2009). Recent advances in assessment methods incorporating spatial structure for fish populations in the Great Lakes have resulted in improved estimates (i.e., reduced bias) and have highlighted the need to reduce exploitation on low productivity stocks to avoid population declines (Berger et al. 2012; Molton et al. 2012, 2013; Li et al. 2014). Movement rates within the Inland Waterway have recently been estimated (Chapter 1), but the management strategy has not yet had time to utilize this information. Therefore, the goal of this paper is to determine the implications of intra-system movement rates on the management of this system, and provide advice to account for this movement. Harvest control rules (or biological reference points) are guidelines that specify acceptable or desirable amounts of catch, fishing effort, or fishing mortality as a function of the current system state (e.g., abundance) estimated from the assessment process (Deroba and Bence 2008). Harvest control rules are designed to help achieve management objectives, and typically incorporate the conservation of a population while allowing for harvest (Deroba and Bence 2008). Constant fishing mortality rate is a harvest control rule that has been widely implemented for the management of freshwater fisheries. For example, many fisheries in the Great Lakes region implement constant fishing mortality rate harvest control rules at a level perceived to pose a low risk of overexploitation (e.g., Staggs et al. 1990; Schueller et al. 2008). However, intermixing, heterogeneous levels of population productivity, and uncertainty (i.e., error) in the assessment and/or implementation of this type of harvest strategy have important management implications (Deroba and Bence 2008; Nieland et al. 2008; Molton et al. 2012, 2013). Uncertainty in the state of fishery systems or their response to management actions is inherent in fishery management decisions (Hayes et al. 2014) and ignoring uncertainty in the management decision process can lead to unsustainable levels of exploitation (e.g., Punt 2006). 28 Specifically for populations assessed using mark-recapture methods, the level of uncertainty for population estimates can be high because the number of fish marked and subsequently examined for marks is often low or the study design lacks the ability to meet model assumptions (Staggs et al. 1990; VanDenAvyle and Hayward 1999). The high level of uncertainty typical of fishery studies and associated risk of overexploitation highlights the importance of explicitly incorporating uncertainty into the management strategy evaluation process (Frederick and Peterman 1995; Katsukawa 2004; Deroba and Bence 2008; Nieland et al. 2008). The purpose of our study was to evaluate the implication of spatial structure on the application of the current harvest policy for walleye populations in a lake-chain where the management strategy does not currently account for intermixing. Our study addresses longstanding concerns of overexploitation on intermixed populations and addresses concerns specific to walleye fisheries. In particular, Rasmussen et al. (2002) highlighted issues associated with setting harvest quotas for walleye fisheries (i.e., spearing and angling) that occur during different time periods that are based solely on spring-spawning population assessments, and therefore overlook post-spawn movements. Our specific objectives were to: 1) evaluate the implications that intermixing had on the allocation of total allowable (TAC) using a constant mortality rate for two fisheries (i.e., spearing and angling) that occur during different time periods, 2) determine the effect that variable levels of spearing mortality and site-specific angling mortality have on the overall exploitation rate for spawning populations that are intermixed during the angling season, and 3) determine the influence that process error for angling season mortality has on spawning population exploitation rates. To accomplish these objectives we took a simulation modeling approached that allowed us to evaluate the risk of overexploitation associated with the use of specific harvest control rules 29 and/or assessment techniques for exploited fisheries (Punt 2006). Simulations have been used to show that incorporating spatial structure results in improved estimates of system state parameters (Porch et al. 2001; Berger et al. 2012; Molton et al. 2013; Li et al. 2014). The strengths and flexibility of using a simulation framework make the approach advantageous for addressing complex issues associated with incorporating spatial structure into management strategy evaluation (Goethel et al. 2011; Molton et al. 2012, 2013; Li et al. 2014). Therefore, we used a stochastic simulation-framework to evaluate the implications of spatial structure on the application of harvest control rules while accounting for plausible and observed levels of exploitation for walleye populations within our study area. Methods Study area Michigan’s Inland Waterway is an interconnected chain of lakes located in the northern Lower Peninsula that consists of four large lakes (Burt, Crooked, Mullett, and Pickerel) interconnected by a series of rivers and smaller tributaries (Figure 2.1). The Cheboygan Lock and Dam on the Cheboygan River, and the Alverno Dam on the Black River located at the northern portion of the Inland Waterway, restrict fish passage and are considered closed to emigration towards Lake Huron or further upstream within the Black River (Figure 2.1). The lakes and rivers of the waterway are oligotrophic, provide various levels of suitable walleye spawning substrate and prey resources, and range from 4.4 km2 (Pickerel Lake) to 70.4 km2 (Burt Lake) in total size (Hanchin et al. 2005a; Hanchin et al. 2005b). Walleye within the Inland Waterway are highly mobile and exhibit spatial structuring that is typical of intermixed populations (Chapter 1). Therefore, we used walleye populations from the Inland Waterway as our model species to evaluate the implications of ignoring movement 30 patterns on the application of harvest control rules for an intermixed population. Recently, Herbst et al. (Chapter 1) estimated among-lake post-spawn movement rates and spawning sitefidelity of walleye populations in the waterway and recommended that the management strategy be evaluated in light of the observed movement patterns. The management objective in the waterway is to maintain self-sustaining walleye populations, while allowing a subsistence spearing harvest during spawning periods and a state-regulated angling harvest from late-April to mid-March. To date there has not been an evaluation of the current management strategy in the waterway addressing the complexities of observed spatial structuring, differences in productivity, and fisheries that occur over varying temporal scales. General approach and model structure We developed a stochastic simulation model to determine the influence of spatial structure on the application of the current harvest control rule for walleye populations in the Inland Waterway. The model was parameterized using lake-specific estimates of abundance (Michigan Department of Natural Resources unpublished data), and estimated post-spawn movement rates, fishing mortality rates, and spawning-site fidelities specific to our study area (Chapter 1; Table 2.1). We evaluated eight scenarios that represented the observed movement and mortality dynamics from walleye populations within the waterway and current management strategies (Table 2.2). These scenarios included various levels and combinations of assessment, implementation, and process uncertainties associated with walleye demographics and harvest management. Specifically, assessment uncertainties for spawning population size at each lake/river and in each year assumed population estimates were unbiased and generated from a normal distribution with known levels of precision (estimated from mark-recapture data). 31 Implementation uncertainty was included for recreational angling mortality in the form of estimated lake-specific fishing mortality rates on summer feeding grounds (Chapter 1) when populations were intermixed. Moreover, total instantaneous mortality rates on summer feeding grounds were allowed to vary in a lake-and-time specific fashion, simulated from a lognormal distribution with a proscribed level of variation (described further below; Chapter 1). We had no direct estimates of tribal spearing harvest, and as such, we conducted a sensitivity analysis to evaluate multiple levels of spearing mortality without including implementation error. For each scenario, we used 1,000 simulations of 200 year projections to ensure that transient dynamics in the initial years did not obscure the distribution of relevant performance metrics. Moreover, because we were interested in the ability of management to meet the maximum 35% exploitation rate for a given year in the presence of movement and uncertainties, we used the value of our performance metrics in the final projection year to construct distributions of harvest performance metrics over all 1,000 simulation iterations. Initial lake-specific abundances were set equal to estimates from a system wide markrecapture survey conducted in 2011 (Michigan Department of Natural Resources, unpublished data). Abundance estimates prior to spearing ( ̂ ) at each location (i) and time (t) were then randomly generated using a truncated normal distribution (to prevent estimates < 0 at low abundances) with the mean equal to the spawning population estimate and the standard errors estimated for each location from the 2011 mark-recapture survey. This allowed us to account for assessment uncertainty when determining the spawning population size for each location, which then affected harvest-policy performance by governing the number of fish being harvested via spearing on the spawning grounds: ̂ 32 . In this equation refers to the population of spawners remaining after spearing, and is the spearing exploitation rate. For our simulation scenarios we used a sensitivity analysis with a set of fixed values for us that ranged from 0.05 to 0.20 to cover levels of spearing exploitation (us) that are currently deemed plausible. The spearing exploitation rates were then held constant across locations for the suite of scenarios evaluated (Table 2.2). The base scenario (i.e., scenario 1) held spearing exploitation constant at us = 0.10, which coincided with the maximum spearing mortality for the current harvest policy in the waterway. Following spearing harvest, the remaining individuals from each lake exhibited postspawn movements to the locations where they experienced recreational angling and natural mortality (Figure 2.2). We used estimated post-spawn movement rates (Chapter 1) to determine the abundance of fish at each location that were available for recreational-angling harvest. Specifically, the matrix of estimated posterior distributions of movement rates from spawning to summer locations were used to generate post-spawn movements ( = movement from spawning site i to summer location j) (Table 2.1). Movement matrices for each of the 1,000 iterations of population projections were randomly drawn from the joint posterior distribution of estimated post-spawn movement rates for walleye in this system (Table 2.1). Thus, the number of fish in a given location j directly after post-spawn movements was simply the sum of all fish moving into that location, where fish moving into the site from a spawning location was calculated as the number of spawners surviving spear fishing multiplied by movement rate from ∑ the spawning to summer location (i.e., ; where = the number of fish at location j directly after post-spawn movements). 33 After post-spawn movements, we removed fish from each lake via natural mortality (M) and lake-specific angling mortality. Specifically, the number of fish surviving angling and natural mortality at a given location j and time t ( ) was where ( ( and ) ) . Where Z is total instantaneous mortality and is equal to the sum of estimated location-specific angling exploitation rates (Chapter 1; summarized in Table 2.1 ) and natural mortality (M=0.3) multiplied by a multiplicative process error ( scenarios that did not incorporate process uncertainty the process errors ( ). For operating model ) equaled zero, whereas scenarios that incorporated process uncertainty in angling and natural mortality allowed for lake- and time-specific variation in mortality during the angling season, with ( ). Multiplicative process errors ( ) were simulated from a normal distribution with a mean of zero and standard deviation that was estimated as the mean of the fitted posterior distribution from previous analyses of walleye movement and mortality in this system ( ; Chapter 1). Our simulation program monitored the spawning population membership at each location directly after post-spawn movement, and thus we could determine total exploitation rates that were spawning-population specific. Process uncertainty was included in three simulation scenarios and was used to evaluate the effects of variation in angling and natural mortality on realized total exploitation rates for each spawning population (Table 2.2). After angling and natural mortality, fish from each spawning population that survived either returned to spawn the next year (t+1) at their location of spawning in the current year (i.e., exhibited spawning-site fidelity), or they remained to join the spawning population of their 34 summer feeding location. Additions to these populations also occurred during population projections in the form of recruitment. The spawning population prior to assessment and spearing harvest at site i at time t+1 was the sum of: 1) survivors over time t that never left site i, 2) fish that spawned at site i at time t but survived at another summer location and then exhibited spawning site fidelity, 3) the fish that spawned in another location at time t but survived the summer at site i and then failed to return to their previous spawning population (and thus joined the population of spawners at site i), and 4) new recruits at site i. The proportions of fish from each spawning population exhibiting spawning-site fidelity ( time, and the vector of ) was assumed constant over values for each of the 1,000 population projections was randomly drawn from the estimated joint posterior distributions of Herbst et al. (Chapter 1). The spawning populations were then projected forward using population-specific fixed recruitment at a level that was consistent with the number harvested from each spawning population during the previous time step, essentially performing a yield-per-recruit analysis (Beverton and Holt 1957). ( ( ) ) This approach for recruitment produced population abundance distributions that provided reasonable approximations for walleye in our study area at the final time step in our simulations, and was consistent with other walleye populations in northern Wisconsin lakes (Schueller et al. 2012). Performance metrics We evaluated study objectives by tracking performance metrics, such as total allowable catch (TAC) values and realized total exploitation rates for each spawning population. These performance metrics allowed us to determine the effect that observed post-spawn movement 35 rates had on the implementation of the current harvest control rule for walleye populations in our study area. We compared scenarios that ignored and scenarios that incorporated observed movement dynamics to assess the implications that movement has on the allocation of harvest under the current harvest control rule for walleye populations within the waterway. The example was used as a performance metric to depict the difference in total allowable catch (TAC) values for the angling fishery by location after accounting for the re-distribution of individuals following the spawning season. We provided location-specific adjustments for the angling season TACs that would be needed to account for post-spawn movement dynamics, while still adhering to the maximum exploitation rate of 35% for each spawning population. Specifically, the angling TAC adjustments were determined by evaluating the difference between scenarios with and without movement, with negative values indicating that fewer fish would need to be harvested to maintain the u = 0.35 for the spawning populations throughout the waterway. Realized total exploitation rates for each spawning population were used as a performance measure to evaluate the implications that movement dynamics and scenarios of differing mortality rates had on the implementation of the harvest control rule (Table 2.2). For each of the scenarios evaluated, realized total exploitation rates for each spawning population was determined by dividing the total harvest from the tribal and angling fisheries at each time t by the spawning population abundance at time t. The total harvest from each spawning population was equal to the sum of spearing and angler harvest. Angler harvest was calculated with the Baranov catch equation (Quinn and Deriso 1999) using estimated angling exploitation rates derived from Chapter 1. Spearing harvest was calculated as the product of the locationspecific spawning population abundance and the spearing exploitation rate (us). For example, the spearing harvest for the Burt Lake spawning population (N=19,464) under a 10% spearing 36 exploitation rate scenario (i.e., us=0.10) is 1,946. Evaluating the differences between the total realized exploitation rates for each of the scenarios allowed us to determine the effect that spatial structure had on the performance and potential risk of overexploitation associated with applying a constant exploitation rate with a maximum u = 0.35 and assuming no movement in the Inland Waterway. Results Evaluations of the base scenarios that represented perfect implementation of the current harvest control rule of (u = 0.35) without (Scenario 1) and with (Scenario 8) process error resulted in annual spawning population exploitation rates that averaged 0.35 at the final time step (t = 200) of the simulation (Figure 2.3). In addition, the location-specific spawning population abundance estimates from the base scenarios were in the range of realistic values based on our initial values and from density estimates of other walleye populations in northern temperate lakes (Table 2.3; Schueller et al. 2012). These realistic and expected outcomes confirmed that our simulation model was operating properly. Asymmetrical post-spawn movement dynamics influenced the abundance of walleye in each location during the post-spawn time period. Therefore to achieve the maximum annual exploitation rate (u = 0.35), the angling total allowable catch (TAC) for each location would need to be adjusted within the waterway (Table 2.3). This result was illustrated by the differences in angling TACs from the base scenario (scenario 1) and a slightly modified version of scenario 1 that excluded post-spawn movement dynamics (Table 2.3). The differences between location-specific angling TACs when movement was included or excluded highlighted how post-spawn movement dynamics influenced the implementation of the harvest control rule 37 in the waterway (Table 2.3). The adjustments in angling TACs needed to implement the maximum exploitation rate and account for the net movement of spawning populations among the locations within the waterway ranged from 280 fish less to 370 more fish (Table 2.3). The greatest adjustments to angling TACs were needed in locations where post-spawn emigration rates or initial abundance estimates were the highest. For example, Mullett Lake’s angling TAC could increase by approximately 370 fish (Table 2.3) because it was a recipient location of individuals from the Burt Lake spawning population that had low post-spawn movement rates but high abundance, and from the Black River spawning population that had high post-spawn movement rates, but low initial spawning abundance (Table 2.1). In contrast, to account for net movement, the angling TACs would need to be decreased in locations that had a net loss (i.e., Burt and Pickerel lakes and the Black River ) of individuals to achieve the maximum angling exploitation rate (ua = 0.25) for each spawning population (Table 2.3). Simulation scenarios that incorporated process error for total instantaneous mortality (Z) (i.e., scenarios 5-8) had no influence on the mean annual exploitation rate for spawning populations under all levels of spearing exploitation and observed fishing mortality rates (Table 2.4). Despite similar mean annual exploitation rates, the variation around those values was much greater when process error was included (Figure 2.3). We used scenarios that included or excluded process error (e.g., scenario 8 and 1, respectively), with all other inputs held consistent, to evaluate the importance of including mortality rate uncertainty. Using this approach we determined that the differences between the mean exploitation rates for individual spawning populations were negligible (Table 2.4). Therefore, we simplified when the analysis of the influence of movement dynamics by focusing on scenarios that excluded process error (i.e., scenarios 1-4) to make inference. 38 Post-spawn movements among the waterbodies in the Inland Waterway influenced rates of exploitation on spawning populations, however, the simulation scenarios representing the current level of estimated fishing mortality and assumed spearing mortality (scenarios 2-4) did not exceed the maximum annual exploitation rate (u = 0.35; Figure 2.4). Given that u = 0.35 is judged to have a low risk of recruitment overfishing, the fishery in the waterway as a whole does not appear to be overfished. The sensitivity analysis that evaluated differing levels of spearing exploitation indicated that total exploitation rates generally did not exceed the maximum from the harvest control rule. We varied spearing exploitation rates (us) from 0.05 to 0.20 during our simulations (scenarios 24), and even at the highest rate of us (i.e., scenario 4) the total exploitation rate (u) for most spawning populations did not exceed 0.35 (Figure 2.4). The exceptions were the Crooked and Pickerel lakes spawning populations, which had mean exploitation rates (u = 0.38 and 0.37, respectively) that were at the upper level of the maximum u when evaluating the scenario with the highest level of spearing exploitation and observed angling u (Figure 2.4). Discussion Implications of spatial structure on fishery implementation Through the use of a simulation model, we illustrated the implications of spatial structure for a constant fishing mortality rate harvest control rule for walleye populations in a lake chain. Asymmetrical post-spawn movement patterns in the waterway influenced the distribution of total allowable catch (TAC) during the angling fishery for each location. In order to achieve maximum exploitation rates (u = 0.35) on individual spawning populations, the TAC values for the angling fishery need to be adjusted to account for intermixing that occurred after the spawning-season assessment of abundance. Similarly, Rasmussen et al. (2002) concluded that 39 when post-spawn movement occurs, the spawning population estimates for individual lakes in a chain are of less use for management of walleye angling harvest that occurs throughout the year in northern Wisconsin lake chains. Despite asymmetrical post-spawn movement patterns and differing angling mortality rates (ua) among locations within the waterway, the estimated total exploitation rate (u) for all spawning populations was below the maximum of u = 0.35. An evaluation of this maximum exploitation rate (u = 0.35) for walleye populations in northern Wisconsin resulted in an extremely low modeled risk of population decline at all values of initial walleye densities (Schueller et al. 2008), which suggests that there is a low risk of overexploitation when assuming that walleye populations from our study system have similar productivity to those evaluated in northern Wisconsin. The work of Schueller et al. (2008), however, was based upon individual lakes and did not account for systems that support intermixed populations or that have variability in productivity. In these types of systems the allowable exploitation rate may differ among locations to minimize the risk for overexploitation. For example, Molton et al. (2012, 2013) evaluated intermixed lake whitefish populations that have variable productivity and determined that low productivity populations were more prone to overexploitation when fished at maximum mortality rates that were deemed sustainable for the higher productivity populations. Therefore, caution should be taken when extending inference from individual lakes that are closed to intermixing and that have consistent productivity because those results may not be applicable for spatially structured populations in lake chains. Observed angling mortality rates within our study systems (Chapter 1) were on lower end of what has been reported from other studies (Schmalz et al. 2011). However, the populations in our study area would approach the maximum exploitation rate (u = 0.35) if angling mortality increased to levels seen elsewhere, and therefore could be at 40 greater risk of overexploitation, especially if higher mortality rates were exerted on low productivity populations. In this study we attempted to incorporate a level of realism by using estimated population demographics (i.e., movement rates, angling mortality, spawning site-fidelity) to simulate population responses to exploitation, but still we were limited in our ability to incorporate recruitment based on empirical estimates of stock productivity. Determining the productive capacity of fish populations is a challenging process, typically requiring relatively long time series of stock and recruitment data. Since such data were unavailable for the Inland Waterway, we projected our populations forward with a fixed number of recruits, essentially performing a yield-per-recruit analysis (Beverton and Holt 1957), instead of using stock-recruitment (S-R) parameters from other walleye populations, which has been the typical approach in other harvest policy evaluation studies. Our decision to use fixed recruitment levels was based upon the desire to avoid the known sensitivity of population demographics to the form of the stock-recruitment relationship (Deroba and Bence 2008). Furthermore, our approach allowed us to address variable levels of productivity among the spawning populations by setting recruitment at a fixed value consistent with the harvest from each population. Therefore, more productive spawning populations that support greater harvest (e.g., Burt Lake spawning population) under the harvest control rule also had a higher level of recruitment. Despite our efforts to address population productivity, however, we recommend developing a stock-recruitment relationship for each of the spawning populations within the waterway to better understand the productive capacity and relative vulnerability of walleye in this system. The addition of realistic S-R parameters would allow managers to build upon our analysis by adding further realism that would benefit future evaluations for other performance measures, such as sustainability, under different maximum 41 exploitation rates (e.g., Nieland et al. 2008; Schueller et al. 2008) or harvest policies (e.g., Deroba and Bence 2008). Addressing sources of uncertainty Uncertainty is innate in fisheries science, with important implications for the evaluation of management strategies (Hayes et al. 2014). Incorporating uncertainty, both in state variables (e.g., abundance) as well as population and fisheries processes, provides more realistic insight into the full response of performance measures to management actions. Most prior work on the effect of uncertainty on harvest policy performance has focused on the influence of errors in population abundance estimates (Deroba and Bence 2008). During this study, we incorporated spawning population abundance uncertainty by simulating the population size each year using the variance from previous assessments. This approach resulted in variability in the abundance point estimate, essentially representing assessment error, which was used to determine harvest allocation. Given the variability in abundance within all simulated scenarios, this study focused less on the effect of abundance uncertainty, but instead focused the effect of process error on the mortality source (Z) during the angling fishery. During these evaluations we found that mean exploitation rates were unchanged under scenarios with and without process error on Z during the angling season. We hypothesize that the influence of process error on a demographic rate, such as Z, has a lesser effect on performance measures than does abundance estimate uncertainty because of the difference in scale. For example, abundance estimates from assessments can be highly uncertain because of challenges in capturing enough fish for precise estimates, as well as potential biases due to violation of model assumptions (Hayes et al. 2007). In contrast, uncertainty associated with estimated demographic rates in this study was less substantial. For 42 example, Herbst et al. (Chapter 1) estimated that the median value of the multiplicative process error for mortality (Z) during the angling season was approximately 0.7 for walleye populations in a lake chain. This estimate of process error resulted in a lesser effect on the population performance measures when compared to the documented greater effects of ignoring population estimate variance (i.e., assessment error) (Frederick and Peterman 1995; Butterworth and Punt 1999; Deroba and Bence 2008; Nieland et al. 2008). Therefore, when evaluating the performance of a constant fishing mortality rate harvest policy, our results suggest the added complexity of incorporating process error on mortality estimates is less imperative than incorporating uncertainty from abundance or biomass assessments. An additional uncertainty in the Inland Waterway revolves around the actual level of exploitation rate during the spearing harvest. The rate of exploitation is currently unknown, however anecdotal information suggests that us of 0.20 or less is a plausible value. Our study addressed this uncertainty by using a sensitivity analysis that explored plausible levels of exploitation that were above and below the current maximum for the spearing fishery. This evaluation indicated that spawning population exploitation rates were less than the maximum mortality rate (u = 0.35 ± 95% CI), even when the maximum spearing exploitation rate was raised to us = 0.20 from the current maximum of us = 0.10. These results, however, only hold true when combined with the current observed level of location-specific fishing mortality rates (ua). Our findings were consistent with other studies reporting that the performance of harvest control rules was not sensitive to the variance around a maximum mortality rate (Sethi et al. 2005; Nieland et al. 2008). All uncertainties, however, add to the realism of describing biological systems and despite performance insensitivities to implementation error, Butterworth and Punt 43 (1999) suggest that all uncertainties have value in evaluations because they can interact with assessment uncertainty and affect performance of harvest policies. Management implications and conclusion The process of managing populations that exhibit seasonal intermixing is individually complex, but the inherent uncertainty of assessing a fishery further exacerbates the difficulties of making harvest management decisions. Two prominent issues associated with managing the spatially structured populations in our study area are: 1) unknown level of spearing harvest implementation and 2) the allocation of angling harvest (ua) on intermixed populations that is set based upon a constant mortality rate (u = 0.35) and an abundance estimate from an assessment that occurs when populations are spatially distinct. This study evaluated these management concerns and determined that under current and plausible conditions of angling and spearing exploitation, each of the spawning populations were being exploited at rates lower than the maximum constant mortality rate. However, to achieve the maximum exploitation rate under perfect implementation of the harvest control rule, the total allowable catch for the angling fishery would need to be adjusted to account for post-spawn intermixing of populations. As such, if these adjustments are not implemented, we recommend that managers focus efforts on monitoring exploitation rates which will ensure that observed exploitation rates remain below maximum levels. Monitoring exploitation rates for desirable sport fish species is critical because these rates are variable (Schmalz et al. 2011), and if they increase management actions will likely be needed to avoid exceeding maximum levels. Furthermore, because the level of productivity influences the sustainability of intermixed populations (Molton et al. 2012, 2013) we recommend that managers develop a stock-recruitment relationship for each spawning 44 population to better understand the variability in productivity among populations in the waterway. If low productivity populations are identified, it may warrant altering maximum exploitation rates to protect them from overexploitation. Finally, based on existing literature and results from this study, it is imperative to include levels of uncertainty for abundance when setting TACs (Butterworth and Punt 1999; Deroba and Bence 2008), but potentially less crucial to include implementation error or process error on demographic rates (This study; Sethi et al. 2005; Nieland et al. 2008). Therefore, we recommend that managers strive to assess abundance during time periods that best represent individual stocks (i.e., spawning periods) and to include the associated level of uncertainty when using those abundance estimates combined with movement rates to evaluate the performance of harvest control rules in the future. In summary, this study evaluated intermixing that re-distributed spawning populations during the angling fishery and determined that total allowable catch (TAC) would need to be adjusted to achieve the maximum constant mortality rate of u = 0.35 to account for this spatial structure. We determined that the mean exploitation rates for individual spawning populations were lower than the maximum exploitation of 35% set by the harvest control rule under the harvest conditions/scenarios we evaluated. Despite these encouraging results, we recommend that managers work towards developing location-specific stock-recruitment (S-R) relationships. An understanding of the S-R relationships would enable managers to evaluate additional performance measures, such as sustainability, without biasing results from the known sensitivities to using stock-recruitment parameters from other populations (Deroba and Bence 2008). Additional value to management would be added by developing a stock-recruitment relationship for each spawning population because it would provide information on productivity, which has implications for sustainable harvest strategies for intermixed populations (Molton et 45 al. 2012, 2013). Lastly, our study provided evidence that the spawning populations were not at risk of overexploitation, if in fact the 35% maximum exploitation rate is associated with a low risk of overharvest in our study area. We do, however, suggest that harvest management proceeds with caution because of the known increases in risk of overexploitation in intermixed fisheries with differing productivity levels. 46 CHAPTER 3 Walleye Foraging Ecology in an Interconnected Chain of Lakes Influenced by Non-Native Species Introduction The introduction of non-native species is a widespread phenomenon that has influenced the ecological processes and interaction of native fish communities (Cucherousset and Olden 2011). Non-native species typically possess life-history traits that allow them to out-compete native species and become prolific within their introduced ecosystem (e.g., Lodge 1993; GarciaBerthou 2007). Established non-native species influence native species at population, community, and ecosystem scales of biological organization (Cucherousset and Olden 2011). Some non-native species, such as dreissenid mussels (Dreissena polymorpha and D. bugensis), have had consistent and well-documented negative ecological effects on native fish communities (Vanderploeg et al. 2002; Pothoven and Madenjian 2008; Campbell et al. 2009; Cucherousset and Olden 2011), while effects from other non-native species (e.g., round goby (Neogobius melanostomus) or rainbow smelt (Osmerus mordax)) have been more variable. The impact of non-native species seems to be context dependent (e.g., Roth et al. 2010; Kornis et al. 2013) and has been linked to the adaptability of native predators to use non-natives for prey (Lumb et al. 2007; Herbst et al. 2013). Few studies, however, have extended evaluations beyond the effect of a single non-native species, particularly for smaller inland lakes (i.e., all lakes smaller than the Great Lakes). This leaves a knowledge gap of how multiple non-native species in the same system will influence the trophic ecology of native predators. Understanding the trophic linkages between non-native and native species will aid in management decisions related to conserving 47 native fisheries, and provide insight into the future productivity of valuable commercial and recreational fisheries. The effect of non-native species on native fish populations might be amplified after they become established and integrated into an ecosystem’s food web. Examples date back to the foundational work of Brooks and Dodson (1965) that documented ecosystem alterations (i.e., trophic cascades) caused by the introduction of alewives (Alosa pseudoharengus) in landlocked lakes. More recently, the diet of a native Great Lakes generalist predator, lake whitefish (Coregonus clupeaformis), has been shown to consume abundant, but less energetically valuable, invasive Dreissena spp., resulting in decreased growth and condition (Pothoven and Nalepa 2006; Lumb et al. 2007; Pothoven and Madenjian 2008). In contrast, zebra mussels were not incorporated into lake whitefish diets in Lake Champlain and thus had little direct influence on their condition or demographics (Herbst et al. 2013), which highlights variability in the response of native to non-native species. In the Great Lakes, studies have historically focused on the influence of individual non-native species, but the focus has recently broadened to evaluate the forage ecology of native piscivores in the presence of trophic shifts and multiple non-native species (Campbell et al. 2009; Kaemingk et al. 2012; Pothoven and Madenjian 2013; He et al. 2014 ). It is unknown, however, whether results from the Great Lakes can be downscaled to inland lakes, which highlights the need to determine how native predators trophic ecology will respond to the presence of multiple non-native species in small inland lakes. The Great Lakes have experienced the introduction and establishment of more than 180 non-native species (Holeck et al. 2004). As such, the Great Lakes are a source for the secondary spread of non-native species to other waters in the region (Vander Zanden and Olden 2008). The increased spread of non-native species from the Great Lakes emphasizes the importance of 48 understanding how these species will influence the trophic dynamics of native fish in smaller inland lakes. As species such as alewife, round goby, and rainbow smelt expand from the Great Lakes, they have the potential to become integrated into the food webs of other waterbodies. If native predators in smaller ecosystems react similarly to round goby as they have in the Great Lakes, round goby have the potential to become a major component in the diet of many predators and thus an important element of the food webs (Dietrich et al. 2006; Madenjian et al. 2011; Kornis et al. 2012; Pothoven and Madenjian 2013). Similarly, the non-native alewife and rainbow smelt became major contributors to predator diets after their introduction into the Great Lakes and other inland lakes and reservoirs (Jones et al. 1994; Krueger and Hrabik 2005; Hobson et al. 2012). Despite the growing knowledge base of how individual non-native species have influenced and become integrated into the food web dynamics of large scale ecosystems (i.e., Great Lakes), few studies have simultaneously evaluated the level of integration of round gobies, alewives, and rainbow smelt into smaller scale food webs. The walleye (Sander vitreus) is a predator that is native throughout the north temperate region, inhabiting many rivers and lakes in North America (Billington et al. 2011). Walleye are a model native predator to determine how pelagic and littoral non-native fish species become integrated into the native food webs of smaller inland lakes. Walleye, like many predators, exhibit ontogenetic dietary shifts from zooplankton to aquatic insects before becoming piscivorous. As piscivores, walleye are generalists preying upon a wide variety of available prey (Chipps and Graeb 2011). Patterns in prey availability and diet of walleye typically follow a seasonal pattern. For example, in the spring when aquatic invertebrates (e.g., Hexagenia spp.) are readily available, these organisms can be a significant portion of the diet (Forney 1974). Walleye typically overlap geographically with yellow perch (Perca flavescens) and as the season 49 continues, walleye utilize this native prey fish as a major dietary item, primarily in late summer and fall when young-of-the-year yellow perch reach a size selected by walleye (Forney 1974, Reeves In review). A major question is how non-native species will affect this prototypical pattern, and the extent to which walleye integrate those non-native species into their diet. In this study, we describe the contribution of several native and non-native forage species to the feeding ecology of walleye within a north temperate lake chain system that has been invaded by multiple non-native species originating from the Great Lakes. Our a priori hypothesis was that walleye diet would be linked to the relative amount of prey in open-water habitats. Two of our study lakes (Mullett and Burt) contain a high percentage of hypolimnetic pelagic habitat, whereas the other two (Crooked and Pickerel) are largely littoral. Therefore, we expected walleye to prey upon the most abundant forage found within the pelagic zone (e.g., yellow perch, alewife, rainbow smelt) in Burt and Mullett lakes, and within the littoral zone (e.g., yellow perch) in Crooked and Pickerel (Chipps and Graeb 2011). Furthermore, we hypothesized that this prey-habitat relationship would be modified in systems with abundant non-native species (i.e., round goby, alewife, and rainbow smelt) potentially drawing predators away from the dominant habitat type. Thus, we expected that the non-native species would be highly integrated into the forage ecology of walleye in Mullett and Burt lakes because these lakes have established populations of round goby and alewife, although alewife have higher densities in Mullett Lake (this study). In contrast, we expected that the diet of walleye in Crooked and Pickerel lakes would be dominated by native littoral prey sources because these lakes have not experienced similar invasions. 50 To evaluate our hypotheses, we used stomach content and stable isotope analyses to (1) quantify the contribution of pelagic versus littoral prey sources for walleye, (2) determine if energy sources differ among lakes based on prey availability, and (3) quantify the contribution of non-native fish species in walleye diets. Stomach content analysis provides direct evidence of species-specific diet on short-term (e.g., hours to days) temporal scales and stable isotopes provide insights on longer-term prey (e.g., multiple months) assimilation patterns (France 1995; Vander Zanden and Rasmussen 1999; Post 2002). Methods Study area The Inland Waterway in the northern portion of Michigan’s Lower Peninsula is a lake chain of four oligotrophic lakes (Burt, Crooked, Mullett, and Pickerel), multiple connecting rivers, and smaller tributaries that drain into Lake Huron (Figure 3.1). The lakes within the waterway range in size, maximum depth, and overall area of habitat types (i.e., littoral and pelagic; Table 3.1). Crooked and Pickerel lakes are similar in size and bathymetry and dominated by littoral habitats (80.3 and 78.8%, respectively; Table 3.1). Burt and Mullett lakes are similar to each other in size and bathymetry, but differ from Crooked and Pickerel in that the pelagic zone makes up the majority (> 68%) of habitat type in these lakes (Table 3.1). Pelagic hypolimnetic habitats in Burt and Mullett lakes support fish assemblages that include cold-water fish species such as brown trout (Salmo trutta) and rainbow trout (Oncorhynchus mykiss), while also supporting cool- and warm-water species in the littoral zone (Hanchin et al. 2005a, Michigan Department of Natural Resources (MDNR) unpublished data). We define the littoral 51 zone as the habitat within the lakes that allows light penetration to the benthos and for our study area is approximately ≤ 4.6 m. The Inland Waterway is connected to Lake Huron via the Cheboygan River (Figure 3.1), which is the suspected pathway for the introduction of alewife, rainbow smelt, round goby, and dreissenid mussels. Abundances of alewife, rainbow smelt, and round goby are unknown, but have all been collected as part of this study and MDNR fish sampling. Despite the paucity of quantitative measures of abundance for round goby, anecdotal angler reports and personal observations provide evidence that this species is abundant in Mullett and Burt lakes, but currently rare/absent in Crooked and Pickerel lakes (MDNR unpublished data). The time frame for the introduction of the pelagic non-native species (rainbow smelt and alewife) is unknown, but was likely in the early 1980s based on anecdotal reports from anglers. Alewife and rainbow smelt seem to be less abundant than round goby throughout the waterway, especially in Crooked and Pickerel lakes based on angler reports and visual observations during this study. Forage assessment Pelagic prey availability was determined using small mesh vertical forage gillnets (FGNs). The FGNs were 3.7m wide and comprised of four panels (0.9m/panel) of variable mesh sizes (9.5, 12.7, 15.9, 19.1mm stretch) and the nets were capable of fishing the entire water column (bottom to surface). These nets were deployed during three sampling events (late Mayearly June, late June-early July, and late July-early August) in 2011 following a depth-stratified random sampling design for each of the four lakes. During each of these sampling events, FGNs were deployed and fished for one net-night in five depth zones that covered all available habitats in each lake (total sets for waterway by depth zone: 0-2.1m = 57, 2.2-6.1m = 95, 6.2-12.2m = 93, 52 12.3-18.3m = 63, > 18.3m = 35). In 2012, FGNs were used during two sampling events (late June-early July, and late July-early August) following the 2011 sampling design. Modifications to the sampling design in 2012 were made to focus efforts on depth zones where sampling was most productive in 2011, so the shallowest depth strata (0-2.1 m) was removed from the sampling design. The selectivity of the FGNs made it difficult to quantify round goby densities because these fish were not vulnerable to the gear. Kornis et al. (2012) evaluated capture methods for round goby and concluded that all methods had limitations. Furthermore, Hayes et al. (2012) describe the complexities with calibrating across gear types. We used FGNs to evaluate prey relative abundance across a variety of species and habitats, but similar to the inefficiencies mentioned in Kornis et al. (2012) we found that round goby were not highly vulnerable to this gear. As such we relied on qualitative observations that were based on personal observation and angler reports. Prey fish collected using the FGNs were identified to species, counted, and measured for total length (mm). Estimates of relative abundance, measured as catch-per-unit-effort (CPUE: fish/net), and species composition of the forage base, except round goby, were calculated for each lake within the Inland Waterway. Stomach content analyses Walleye diets were sampled using gastric lavage on individuals incidentally collected during forage gillnetting and fall juvenile electrofishing assessment sampling and from whole stomachs of harvested walleye donated by anglers during 2011-2013. Stomach contents collected with gastric lavage and whole stomach samples from anglers were individually labeled and stored frozen. In the laboratory, stomachs were dissected and all contents were removed, identified to species when possible, counted, and weighed (g). All contents were then classified 53 into the following dominant prey categories: aquatic insects, yellow perch, round goby, alewife and rainbow smelt, crayfish, unknown fish, and other. The observed diet was quantified as percent composition by wet weight (g) for each of the categories by season (spring, summer, fall/winter) and lake (Mullett, Burt, Crooked, and Pickerel). Empty stomachs were included in the analysis and counted as zeroes. The spring season was defined as April to June, the summer season as July to September, and the fall/winter season as October to March. Stable Isotope analyses Stable isotope analysis is beneficial for determining long-term prey assimilation patterns because δ13C values tend to be conserved from prey to predator, thus providing insight on energy source, which can be derived from littoral (e.g., attached algae) or pelagic (e.g., phytoplankton) production. The ability to decipher energy sources is associated with the enrichment of the δ13C of the base of the littoral food web relative to the base of the pelagic food web (France 1995; Vander Zanden and Rasmussen 1999; Post 2002 ). Muscle tissue samples from multiple individuals of various species were collected from each of the lakes within the study area from late-July to mid-August, 2012 (Appendix 3.1). Dorsal muscle tissue was collected and used for analysis from fish species large enough in size to attain a fillet (i.e., walleye and yellow Perch). For all other taxa, we homogenized whole body samples after removing the head, stomach, and digestive tract. Muscle tissue samples were used because they provide a long-term (i.e., multiple months) integrated image of the isotope composition of the food consumed during the growth period (Perga and Gerdeaux 2004). A wide size range of walleye (335 – 600 mm TL) were collected, whereas only yellow perch ≤ 150 mm TL (range = 82 - 150 mm) were analyzed because that size range represents the size range 54 that is most susceptible to predation from walleye (Overman and Parrish 2001; Reeves in review). The isotopic values of a consumer alone is not generally sufficient to infer trophic position or carbon source (Post 2002), so we collected zooplankton and gastropods to represent the isotopic baselines from the pelagic and littoral food webs, respectively. We quantified stable isotope values for gastropods, zooplankton, and Hexagenia spp. using composite samples comprised of 5 to 25 individual samples from each taxa (Appendix 3.1). For example, an individual gastropod stable isotope sample consisted of 25 homogenized individuals. Individual or composite samples were dried, a subsample of approximately 1.25 mg dry of each sample was placed in 3 x 5 mm tin capsule, labeled and analyzed for dual δ13C and δ15N analyses at the Stable Isotope Facility at the University of Wyoming using a Costech 4010 Elemental analyzer coupled with a Thermo Delta Plus XP IRMS. Stable isotope values are conveyed in δ notation where δ13C or δ15N = [(Rsample/Rstandard)-1] x1000, where R is 13C/12C or 15N/14N. Carbon and nitrogen were standardized using Pee Dee Belemnite carbonate and atmospheric nitrogen, respectively. Precision was ± 0.05 for nitrogen and ± 0.1 for carbon based on laboratory standards. We quantified trophic position for each species analyzed using the equation from Vander Zanden and Rasmussen (2001). ( ) where 3.4 is the assumed trophic fractionation constant between trophic levels (Minagawa and Wada 1984) and the value 2 is added because trophic position was estimated relative to primary consumers rather than producers (Post 2002). Within the equation to calculate trophic position, we used walleye as the δ15Nconsumer and δ15Nbaseline was zooplankton, which for this study was a composite sample of non-predatory zooplankton species (i.e., seston). The contribution of 55 pelagic and littoral prey sources to walleye energy intake was determined using a δ13C mixing model (Vander Zanden and Rasmussen 2001). where the subscript C indicates the consumer, P is the pelagic end member, and L is the littoral (benthic) end member. The specific end members were chosen because zooplankton represent the primary producers and form the base of the pelagic food web and gastropods represent the primary producer and base of the littoral food web (Post 2002). These species were used for calculating trophic position and percent littoral, respectively, for individual species within all four lakes of the waterway. Using these two organisms as baselines for the two mixing models constrained estimates of trophic position to be ≥ 2.0 and percent littoral to be between 0 and 100%. We found that the majority of the lake specific isotope data adhered to the assumption of normality (Shapiro-Wilk Normality test); therefore, we used one-way ANOVA to evaluate differences in isotopic signatures (δ13C and δ15N) of the pelagic and littoral baseline species (zooplankton and gastropods, respectively) among the four lakes within the waterway. We found statistical differences in δ13C and δ15N values of our baseline species between lakes, and therefore we used percent littoral and trophic position instead of raw δ13C and δ15N values. After standardizing isotopic values, we used a general linear model (i.e., ANCOVA) to determine if walleye trophic position was influenced by total length for each location, which has been reported for walleye (Overman and Parish 2001). The isotopic niche size takes into account the nitrogen range and carbon range, providing information on the breadth of energy sources being consumed, as determined by the δ13C range and also provides insight into the trophic level feeding patterns in the food web by describing the 56 δ15N range utilized by a population (Layman et al. 2007; Jackson et al. 2011, 2012). We tested for differences in the isotopic niche size of walleye among using MANOVA. To further assess differences among walleye isotopic niches, we used the Stable Isotope Bayesian Ellipses in R (SIBER) package within R (R Development Core Team 2012). Isotopic values from walleye sampled in each of the lakes were used for these calculations. Detailed methodology and calculations for population metrics described by Jackson et al. (2012) were followed to calculate standard ellipse areas corrected for small sample sizes (SEAc). The SEAc values were then used to evaluate the isotopic niche size among walleye populations within the waterway. Because we compared isotopic niche among walleye populations, we used estimates of the trophic position and the percent littoral for individual walleye to compare the isotopic niche for walleye between lakes. We used a non-parametric jack-knife resampling technique to determine SEAc differences among lakes instead of using the Bayesian bootstrap methodology described by Jackson et al. (2011). Lastly, we used the jack-knife SEAc estimates and 95% confidence intervals to evaluate differences in isotopic niche (i.e., SEAc) of walleye within the waterway. Results Forage assessment Yellow perch was the dominant prey species collected using FGNs, consisting of approximately 85% of the total catch (Table 3.2). The overall CPUE differed among lakes and varied by year (Table 3.2); differences in overall CPUE from 2011 to 2012 were primarily linked to the change in yellow perch catch rates. For example, the overall CPUE in Mullett Lake decreased from 2011 to 2012, but most of that change is due to the decrease in yellow perch CPUE (Table 3.2). The primary non-native species collected in the FGNs were alewife, and to a 57 lesser extent rainbow smelt and round goby. The FGN nets were not effective for sampling round goby, however, during stable isotope sampling efforts we observed high angling capture rates (> ~35/hr; personal observation) and these fish were visibly abundant in nearshore rocky substrate within Mullett and Burt lakes. Mullett Lake had the highest two year average catch rate (12.9 fish/net) and also had the highest catch rates of non-native fish (Table 3.2). Alewife relative abundance in Mullett Lake increased during the second year of the study and was the second most abundant species in the catch (Table 3.2). Pickerel Lake, the furthest from Lake Huron, was the only lake where non-native species were not collected as part of this study (Table 3.2) or reported through other means (e.g., angler reports or agency sampling). Stomach content analysis A total of 1,484 walleye stomachs were collected to quantify diet within the Inland Waterway (Table 3.3). Throughout the three-year collection period, the sample size of stomachs turned in from anglers increased from all locations and seasons and allowed for seasonally (spring, summer, fall/winter) stratified diet analysis. Despite our outreach efforts, there were still low sample sizes that limited the precision of estimated diet from Crooked, Mullett, and Pickerel lakes in the fall/winter season (Table 3.3). Walleye exhibited seasonal foraging patterns that integrated non-native prey fish into their diets. Walleye fed primarily on aquatic insects, in particular Hexagenia spp., in the spring and transitioned into a fish or crayfish dominated diet in the summer and fall/winter (Table 3.3). The exception to this seasonal pattern was walleye from Mullett Lake where yellow perch dominated the diet during the spring (Table 3.3). The diet of walleye from Crooked and Pickerel Lakes contained fewer taxa than Mullett and Burt lakes and consisted of primarily crayfish (> 58 50% in the summer) with other prey fish (i.e., miscellaneous Notropis spp.) and insects composing the remainder of the diet from walleye in those lakes (Table 3.3). Non-native prey species, where present, were observed in walleye diets. For example, round goby comprised the largest portion of the diet during the summer and fall seasons in Mullett and Burt lakes (Table 3.3). In Mullett Lake, for example, round goby made up approximately 77% of walleye diet during the fall/winter season, although low sample size limits the precision of that value. Similarly, round goby in Burt Lake were heavily utilized in the summer and fall/winter, making up 36.4% and 42.9% of the diet, respectively, during those seasons (Table 3.3). Pelagic non-native prey species (i.e., alewife and rainbow smelt) were preyed upon to a lesser extent by walleye than we anticipated (Table 3.3). In fact, only walleye from Mullett Lake contained alewife and/or rainbow smelt in their stomachs, and those two species represented only a low percentage of the diet in the summer (6.9%) and fall/winter season (2.9%) within Mullett Lake. Crayfish were the primary diet item in Crooked and Pickerel lakes in the summer season, contributing approximately 50% of the total weight of prey items (Table 3.3). Although most crayfish in the stomachs were difficult to decipher to species because of digestion, we confirmed the presence of non-native rusty crayfish (Orconectes rusticus) in walleye from Crooked and Pickerel lakes. Stable isotopes The isotope values (δ13C and δ15N) of the pelagic (zooplankton) and littoral (gastropods) baseline organisms were significantly different among lakes (pelagic baseline ANOVA: F3,15 = 15.93, P < 0.001; littoral baseline ANOVA: F3,15 = 9.96, P < 0.001). Therefore, comparisons of isotopic niche of walleye within the waterway were evaluated using the relative measures of 59 trophic position and percent littoral. The trophic position of walleye was significantly and positively related to length, but there was a significant length by location interaction term (F3,70= 3.93; P = 0.0118). As such, we evaluated the implications of using a site-specific regression to adjust lake-specific measures of trophic position to a common length (412 mm). Our evaluation indicated that there was a biologically negligible effect on the outcome and as such, we focused our interpretation of isotopic analyses on observed (i.e., unadjusted) mean values. (Appendix 3.2). The lake-specific community stable isotope values illustrated the likely prey items of walleye within the waterway. For example, the percent littoral value of walleye within Burt Lake was similar to that of the round goby (Figure 3.2), which was also consistent with diet data that indicated this non-native was an important source of energy (Table 3.3). Similarly, the percent littoral value for walleye within Pickerel and Crooked lakes was similar to crayfish, which corresponded with the primary item observed in their diet (Figure 3.2). Burt and Mullett Lake walleye exhibited a different isotopic pattern, apparently relying on a variety of prey items that range in their energy source (i.e., percent littoral value), as shown by a percent littoral value that envelopes a range of prey items, including round goby, yellow perch, and crayfish (Figure 3.2, Figure 3.3). Mullett Lake walleye, however, had an average percent littoral value greater than Burt Lake which indicated these fish relied more on littoral prey (i.e., round goby and crayfish). The isotopic niche size and energy source of walleye differed significantly among lakes (MANOVA: F3,74 = 22.36 , P < 0.001) and illustrated differences in lake-specific variability among individuals. The SEAc and associated confidence intervals indicated that Mullett Lake walleye had the largest isotopic niche size, followed in decreasing order by Pickerel Lake, Burt Lake, and Crooked Lake (Table 3.4). The SEAc results suggest that the large niche size for 60 Mullett Lake walleye was related to feeding on prey from a wide range of energy sources (Figure 3.3; Table 3.4). Isotopic values of individuals from Pickerel Lake exhibited similar patterns of variability in energy sources, ranging from primarily pelagic to primarily littoral (Figure 3.3). The variability in trophic position associated with Pickerel Lake individuals likely explains the increased estimate of SEAc (Table 3.4). In contrast, walleye from Burt and Crooked lakes had a small niche size, which suggests that individuals were attaining energy from a less diverse assortment of prey. The majority (< 95%) of individual walleye from Crooked Lake, for example, relied exclusively on littoral prey (Figure 3.3, Table 3.4). Discussion In this study, we show that a native predator with a generalist feeding strategy readily incorporated non-native prey items into their trophic ecology in small lakes, which was consistent with the predatory responses observed in larger systems, such as the Great Lakes (Dietrich et al. 2006; Madenjian et al. 2011; Kornis et al. 2012; Pothoven and Madenjian 2013). Although the establishment of non-native species can alter prey communities and drastically influence historical food web structure and energy pathways (e.g., Campbell et al. 2009; Roseman et al. 2014), the ability of walleye to shift their feeding habits to include non-native species is one way the impacts of non-native fishes can be reduced. Prey community Native and non-native fish species co-occurred in a subset of our study lakes, providing a diverse prey base for native predators. Although we were unable to quantify the density of round goby in the study area, angler reports and visual observations (this study) indicated that this non- 61 native benthivore was abundant within Mullett and Burt lakes, but not apparent within Crooked and Pickerel lakes. Although round goby were present within two of our study lakes, there was no evidence suggesting round goby establishment has led to a decline in native fish species as has been observed in the Great Lakes (Kornis et al. 2012). Furthermore, despite establishments of alewife and rainbow smelt in Mullett Lake and alewife in Burt Lake, the native yellow perch was the most abundant fish species collected during our study. This finding was in contrast to other systems where the decline of native species abundance has been linked to the establishment of the non-native alewife and rainbow smelt (Beisner et al. 2003; Mercado-Silva et al. 2007; Madenjian et al. 2013). The significance of the co-occurrence of native and non-native forage species (i.e., alewife, rainbow smelt, round goby, and yellow perch) is that it provides multiple prey resources for native predators in the littoral and pelagic habitats within our study lakes. Foraging Ecology Stomach content analysis Based on previous studies, we expected walleye to rely on yellow perch as a primary prey (Chipps and Graeb 2011), but our results in the Inland Waterway were not consistent with this expectation. Yellow perch were preyed upon less than would be expected, based on their high relative abundance in our study lakes and the generalist feeding strategy that is typical for walleye. Instead, our observations were similar to the predatory responses reported in the Great Lakes, where round goby are abundant and integrated into the diets of native predators. Specifically, our results were similar to Reyjol et al. (2010), who documented that round goby had a higher contribution to predator diets than yellow perch, which were the most abundant native prey fish in the St. Lawrence River system. Furthermore, walleye unexpectedly relied on 62 crayfish as a dominant prey source within Crooked and Pickerel Lakes. To date, few studies have documented such a high reliance on crayfish as prey items (Chipps and Graeb 2011); therefore, the integration of crayfish into the feeding patterns of walleye in our study lakes was unexpected and extends the knowledge of foraging strategies for this predator. Interestingly, Crooked and Pickerel lakes had the lowest forage fish CPUE, suggesting that the reliance on crayfish might have been related to forage fish limitations. Consistency between stomach content and stable isotope analysis We did not observe a consistent match between isotopic values and observed diets in all of our study lakes. Specifically, isotopic signatures and observed diets were inconsistent for Crooked and Pickerel lakes, but were similar for Burt and Mullett lakes. In Crooked Lake the observed diet varied by season with the dominant prey resources originating from a breadth of littoral (i.e., crayfish) and pelagic (i.e., Hexagenia spp. and fish species) energy sources. Typically a broad range of dietary items results in a larger isotopic niche size (e.g., Lodge 1993; Guzzo et al. 2011), but this pattern was not observed for walleye in Crooked Lake. Similarly, the isotopic niche size disagreed with expectations based on observed diet of walleye in Pickerel Lake. The observed diet of walleye in Pickerel Lake was narrow, with crayfish as the dominant diet item during spring and summer – the only two seasons with available data. This single taxa focused feeding strategy led us to anticipate a small isotopic niche size, which was not observed. The inconsistencies between the results from the two methods likely stems from the difference in their effective temporal scales. Additionally, differences between observed diet and isotopic signatures could be the result of low seasonal sample sizes of observed diets from locations within the waterway. 63 Relationship between foraging strategy and habitat The foraging ecology of walleye generally did not match the dominant habitat type available, especially when abundant non-native prey (i.e., round goby) inhabited the less prevalent habitat types. This finding agreed with our hypothesis that the presence of abundant non-native species (i.e., round goby) would influence habitat choice in each lake. In the two lakes with extensive littoral zones (Crooked and Pickerel) and low density of non-native species, the diet of walleye indicated they primarily chose littoral habitats. In contrast, Burt and Mullett Lakes have a large pelagic zone supporting populations of pelagic prey fishes such as alewife, rainbow smelt, and yellow perch. Although yellow perch were heavily consumed in Mullett Lake and to a lesser extent in Burt Lake, alewife and rainbow smelt were not frequently utilized by the native predator, despite alewife being the second most abundant prey species collected during the forage assessment in Mullett Lake. Similarly in Burt Lake, walleye did not heavily depend on pelagic energy, but instead utilized a greater portion of littoral prey (i.e., round goby). This result is in contrast to findings from other systems with substantial pelagic habitats and prey resources. In other lakes supporting alewife and rainbow smelt, these two species typically account for a large portion of walleye diets (Jones et al. 1994; Overman and Parrish 2001; Krueger and Hrabik 2005; Fielder and Thomas 2006; Hobson et al. 2012). Furthermore, our results were contrary to Fielder and Thomas (2006) who reported that alewife were more common in walleye diets from Saginaw Bay, a location where alewife and yellow perch co-occurred. In Mullett and Burt lakes, we suspect that the abundance and vulnerability of round goby was great enough to draw walleyes away from the extensive pelagic zone of these lakes, explaining, at least in part, the small contribution of alewife to the diet of walleye in these lakes, and the limited reliance of prey from the dominant habitat type. 64 The establishment of non-native prey fish species within Mullett and Burt lakes provided additional prey sources; however, this increase in available prey types was not consistently reflected in the patterns of trophic niche size. Although the trophic niche size of Mullett Lake walleye was the largest, walleye from a lake with similar habitat and invasion history (i.e., Burt Lake) was substantially lower. On the other hand, Pickerel Lake, a lake without established populations of pelagic or littoral non-native fish prey species, had a large isotopic niche size. We hypothesize that the trophic niche size is a function of at least two major factors. The availability and utilization of diverse prey and habitat types within a lake is one factor that would be expected to increase niche size by broadening the base from which fish draw their energy resources. The second factor is the degree to which individual fish specialize in either prey types or habitats. Specialization among individuals would increase the variance in trophic utilization among fish, broadening the niche size of the population as a whole. The large isotopic niche size and breath of percent littoral and trophic position values of walleyes from Mullett and Pickerel lakes suggest to us that the degree of specialization among individual fish may vary across systems (e.g., Warburton et al. 1998; Herbst et al. 2013). These results were unexpected and suggest that measures of niche width may provide insight into habitat choice and individual variation in how predators integrate pelagic and littoral prey into their foraging ecology. Management implications of non-natives in the food web The implications of the inclusion of non-native prey in native predator diets in the Inland Waterway have yet to be determined. The inclusion of round goby into walleye feeding has been associated with increased growth rates for some predators in the Great Lakes (Steinhart et al. 2004; Pothoven and Madenjian 2013). The proposed mechanism for increased growth pattern is 65 that round goby act as a conduit of previously unusable prey resources. Round goby consume large quantities of dreissenid mussels and therefore have the ability to convert energy that would otherwise be sequestered in dreissenid biomass (Johnson et al. 2005). Round goby consumption of dreissenid mussels is significant because smaller inland lakes within North America have experienced introduction and establishment of zebra mussels (Benson 2014). Thus, round goby have the potential to provide energy sources to higher trophic levels in dreissenid-invaded systems, thereby potentially offsetting some of the negative effects of these non-native mussels. Although evaluating growth dynamics was outside the scope of this study, unpublished data from the Michigan Department of Natural Resources indicates that walleye in Crooked and Pickerel lakes grow more slowly than in Burt and Mullett lakes. Dreissenid mussels are established in Crooked and Pickerel lakes, and it seems likely that the introduction of round goby may increase growth rate of walleye. The diet and isotopic values of Crooked and Pickerel lakes walleyes reveal that crayfish are the primary littoral energy source being utilized. However, round goby have a higher energy density (4,240 J/g in Pothoven and Madenjian 2013) than do crayfish (3,766 J/g in Roell and Orth 1993), and we hypothesize that if round goby become established in Crooked and Pickerel lakes, walleye diets will shift from crayfish to round goby, and their growth rates will increase assuming consumption rates remain the same. Increased growth rates associated with feeding on round goby, however, have not been consistent for all species or across all locations evaluated (Vanderploeg et al. 2002; Pothoven and Madenjian et al. 2013). Although positive effects on growth have been observed in some situations (Jones et al. 1994; Steinhart et al. 2004; Pothoven and Madenjian 2013), managers need to also consider the potential deleterious effects on native fish population demographics. Adverse effects of nonnative species have been variable in extent and magnitude within the Great Lakes, but highlight 66 potential problems. For example, while alewives have provided a suitable prey resource for adult piscivores in the Great Lakes, this species has limited the reproductive success of many native fish populations throughout the Great Lakes (e.g., Madenjian et al. 2013). In addition, round goby are known egg predators in the Great Lakes and therefore at high densities have the ability to negatively impact the recruitment of native fish populations (Chotkowski and Marsden 1999; Steinhart et al. 2004). Results from available inland studies, however, suggest that the magnitude of the effect on the native predator’s population demographics seems to be related to the predator’s adult stock biomass and ability to exert a significant level of predatory control on the non-native pelagic forage base (Krueger and Hrabik 2005; Roth et al. 2010). For example, after becoming established in many northern Wisconsin lakes, rainbow smelt have had a negative influence on the recruitment dynamics of walleye (Mercado-Silva et al. 2007), despite being heavily consumed by adults (Krueger and Hrabik 2005; Roth et al. 2010). Understanding these types of population responses is imperative because non-natives can directly influence the abundance of native prey and therefore have the potential to affect the dynamics of native predators dependent upon them. Conclusion Native predators in smaller inland lakes exhibit flexibility in their response to non-native prey species. Similar to predators in the Great Lakes, walleye from our study lakes have integrated round goby, where present, into their forage ecology. Our results also indicated that despite having accessible non-native pelagic forage (i.e., alewife and rainbow smelt in Mullett and alewife in Burt), walleye had limited usage of that prey source. Consequently, our hypothesis that foraging ecology would be linked to dominant habitat type was not supported. 67 Furthermore, substantial variability in the foraging ecology of walleye across our study lakes highlights the generalist feeding strategy common for walleye. Our approach and findings indicate that introduced and established non-native prey species, round goby in particular, provide an energy source for native predators and provides insight to the potential response of native predators in other small inland lakes. Furthermore, these findings allow resource managers to have a better understanding of how native fish populations will respond to and integrate nonnative species into the native food web as the secondary spread of these invaders from the Great Lakes continues. Whether newly introduced non-native species serve as prey, or remain unutilized, depends on the degree to which native predators can adapt to incorporate these new potential prey sources into their diets (e.g., Fuiman and Magurran 1994; Strakosh and Krueger 2005). We suggest that predator response to new species introductions is likely to be context dependent, and warrants further investigation for multiple systems and predators to determine the full extent of how non-native species are integrated into the food web and their influence on native communities. 68 APPENDICES 69 APPENDIX 1.0 Supporting tables and figures for Chapter 1: A State-Space Model for Estimating Walleye (Sander vitreus) Movement and Fishing Mortality in Interconnected Lakes. Table 1.1: Model set representing the multiple hypotheses evaluated to represent post-spawning walleye movement and demographics in the Inland Waterway, 2011-2013. Combinations of lake specific and time varying parameters for movement (φ), spawning-site fidelity (ψ), and fishing mortality (F) were evaluated using Deviance Information Criteria (DIC). F(.) is constant fishing mortality for each lake and time. The best fit model was also modified and fit without process error (*) to evaluate model support. Model Number 1 3 4 2 6 8 5 7 1* Structure φ(lake), ψ(lake), F(lake) φ(lake), ψ(lake), F(lake*time) φ(lake), ψ(lake), F(.) φ(lake), ψ(lake), F(time) φ(lake*time), ψ(lake), F(lake) φ(lake*time), ψ(lake), F(lake*time) φ(lake*time), ψ(lake), F(.) φ(lake*time), ψ(lake), F(time) φ(lake), ψ(lake), F(lake) without process error DIC 461.3 467.8 488.6 491.1 516.6 522.0 536.2 538.3 573.3 70 Delta DIC 0.0 6.5 27.3 29.8 55.3 60.7 74.9 77.0 112.0 Table 1.2: Location specific post-spawning movement rates (with 95% credible intervals) estimated by models 1 and 3 under the base assumption for tag shedding rate (Ω = 0.14). Spawning Location Black River Mullett Lake Burt Lake Crooked Lake Pickerel Lake Black River 0.18 (0.04, 0.48) 0.07 (0.0, 0.29) 0.00 (0.0, 0.01) 0.01 (0.0, 0.04) 0.01 (0.0, 0.06) Spawning Location Black River Mullett Lake Burt Lake Crooked Lake Pickerel Lake Black River 0.14 (0.04, 0.45) 0.06 (0.01, 0.34) 0.00 (0.0, 0.01) 0.01 (0.0. 0.03) 0.01 (0.0, 0.06) Model 1 Feeding Location Mullett Lake Burt Lake 0.77 (0.48, 0.92) 0.03 (0.0, 0.08) 0.86 (0.63, 0.94) 0.06 (0.02, 0.12) 0.08 (0.05, 0.10) 0.90 (0.87, 0.93) 0.00 (0.0, 0.02) 0.06 (0.03, 0.09) 0.01 (0.0, 0.04) 0.09 (0.05, 0.14) Model 3 Feeding Location Mullett Lake Burt Lake 0.81 (0.50, 0.92) 0.03 (0.0, 0.09) 0.86 (0.58, 0.94) 0.06 (0.02, 0.12) 0.07 (0.05, 0.10) 0.91 (0.88, 0.93) 0.00 (0.0, 0.02) 0.05 (0.03, 0.09) 0.01 (0.0, 0.04) 0.08 (0.04, 0.14) 71 Crooked Lake 0.01 (0.0, 0.03) 0.01 (0.0, 0.02) 0.01 (0.0, 0.01) 0.87 (0.81, 0.91) 0.19 (0.13, 0.25) Pickerel Lake 0.01 (0.0, 0.05) 0.01 (0.0, 0.03) 0.01 (0.0, 0.01) 0.06 (0.03, 0.10) 0.70 (0.61, 0.80) Crooked Lake 0.01 (0.0, 0.03) 0.01 (0.0, 0.02) 0.02 (0.01, 0.02) 0.87 (0.82, 0.92) 0.19 (0.13, 0.25) Pickerel Lake 0.01 (0.0, 0.05) 0.01 (0.0, 0.03) 0.01 (0.0, 0.01) 0.06 (0.03, 0.10) 0.70 (0.62, 0.78) Table 1.3: Fishing mortality rates (with 95% credible intervals) estimated from model 1 with and without process error using the base assumption of an instantaneous tag shedding rate of 0.14. Sensitivity of fishing mortality rates estimated using model 1 (best fit model) with three different instantaneous tag shedding rates (Ω = 0.04, 0.14, and 0.24). Fishing mortality rates Location Black River Mullett Lake Burt Lake Crooked Lake Pickerel Lake Model 1: w/ process error Ω = 0.04 Ω = 0.14 Ω = 0.24 0.15 (0.02, 0.29) 0.15 (0.02, 0.28) 0.15 (0.02, 0.29) 0.18 (0.13, 0.26) 0.18 (0.13, 0.26) 0.18 (0.13, 0.26) 0.19 (0.16, 0.23) 0.19 (0.16, 0.23) 0.19 (0.16, 0.23) 0.31 (0.26, 0.37) 0.31 (0.26, 0.37) 0.32 (0.26, 0.37) 0.28 (0.21, 0.36) 0.28 (0.21, 0.36) 0.28 (0.21, .036) 72 Model 1: w/o process error Ω = 0.04 Ω = 0.14 Ω = 0.24 0.15 (0.02, 0.28) 0.14 (0.02, 0.28) 0.13 (0.02, 0.28) 0.19 (0.14, 0.26) 0.18 (0.13, 0.26) 0.17 (0.13, 0.26) 0.21 (0.17, 0.25) 0.20 (0.17, 0.23) 0.19 (0.17, 0.22) 0.37 (0.30, 0.43) 0.34 (0.28, 0.40) 0.33 (0.28, 0.38) 0.25 (0.19, 0.28) 0.23 (0.18, 0.30) 0.23 (0.17, 0.29) Table 1.4: Location specific post-spawning movement rates estimated from model 1 with and without process error using different assumed tag shedding rates (Ω = 0.04, 0.14, and 0.24). Black River Spawning Location Black River Mullett Lake Burt Lake Crooked Lake Pickerel Lake 0.04 0.19 0.08 0.00 0.01 0.01 Spawning Location Black River Mullett Lake Burt Lake Crooked Lake Pickerel Lake Black River 0.04 0.14 0.24 0.17 0.19 0.19 0.06 0.07 0.07 0.00 0.00 0.00 0.01 0.01 0.01 0.01 0.01 0.01 0.14 0.18 0.07 0.00 0.01 0.01 0.24 0.18 0.06 0.00 0.01 0.01 Mullett Lake 0.04 0.77 0.85 0.08 0.00 0.01 0.14 0.77 0.86 0.08 0.00 0.01 Model 1 Feeding Location Burt Lake 0.24 0.78 0.86 0.08 0.00 0.01 Crooked Lake 0.04 0.14 0.24 0.04 0.14 0.24 0.03 0.03 0.03 0.01 0.01 0.01 0.06 0.06 0.06 0.01 0.01 0.01 0.90 0.90 0.90 0.01 0.01 0.01 0.06 0.06 0.06 0.87 0.87 0.87 0.09 0.09 0.09 0.19 0.19 0.19 Model 1: w/o process error Mullett Lake Burt Lake Crooked Lake 0.04 0.14 0.24 0.04 0.14 0.24 0.04 0.14 0.24 0.78 0.77 0.77 0.02 0.02 0.02 0.01 0.01 0.01 0.87 0.86 0.86 0.06 0.06 0.06 0.01 0.01 0.01 0.08 0.08 0.08 0.89 0.89 0.89 0.02 0.02 0.02 0.00 0.00 0.00 0.05 0.05 0.05 0.88 0.88 0.88 0.02 0.02 0.02 0.09 0.09 0.09 0.19 0.19 0.19 73 Pickerel Lake 0.04 0.01 0.01 0.01 0.06 0.70 0.14 0.01 0.01 0.01 0.06 0.70 0.24 0.01 0.01 0.01 0.06 0.70 Pickerel Lake 0.04 0.14 0.24 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.06 0.06 0.06 0.69 0.69 0.69 Figure 1.1: Map of northern Michigan’s Inland Waterway that consists of four lakes (Burt, Crooked, Mullett, and Pickerel) and four major connecting rivers (north to south through the lakes: Cheboygan River, Black River, Indian River, and Crooked River). 74 Figure 1.2: Conceptual model depicting the process how a single cohort (e.g., Burt Lake cohort 1) is tracked through time using the presented tag-recovery model. For example, after the initial tagging, which coincides with the spawning period, each individual within Burt cohort 1 has the ability to move to any location within the waterway or can remain in Burt Lake. Following that post-spawn movement the individuals then experience the population and observation processes that are representative of the location they moved to after spawning. Prior to time step t+1, individuals either exhibit spawning-site fidelity and return to their original tagging location (i.e., Burt Lake) or remain in the location they emigrated to. Following the spawning period those individuals once again have the ability to move freely throughout the waterway. 75 Figure 1.3: Comparison of observed tag recoveries (red line) and the predicted distribution of tag recoveries with Bayesian p-values for each cohort during 2011-2013. 76 Figure 1.4: Posterior (and prior) distributions of location specific fishing mortality rates (F) obtained from the best fit model with the base assumption for instantaneous tag shedding rate (Ω = 0.14). 77 APPENDIX 1.1 The explanation for the derivation of prior distribution for fishing mortality ( ̂ ). Pooled catch curve analyses provided estimates of ̂ and ̂ . The estimate of Z is an approximately normally distributed random variable; thus instantaneous fishing mortality is a linear function of a normal random variable ( ̂ Rice (2007; pg. 59): If ( ) and of instantaneous fishing mortality we assumed ̂ ( , then ; thus ( ̂ ) From ). To derive a common prior distribution for estimates and ). 78 , and therefore ̂ (̂ )→ APPENDIX 1.2 Number of individuals released and recovered by location and year for each cohort. Table 1.5: Number marked and recovered by location and year for the 15 tag cohorts, which was used to inform the tag-recovery model. Cohorts 1-3 for each location correspond to the individuals tagged during spring-spawning in 2011-2013, respectively. Recovery years correspond to the annual fishing season that the fish were captured and returned. For example, cohort 1 from Burt Lake was tagged during spring-spawning in 2011 and the recovery years 1-3 correspond with the number of tagged individuals recovered and reported during the 2011-2013 fishing seasons. cohort 1 (n= 5,468) Recovery location Burt Mullett Crooked Pickerel Black River Recovery year 1 2 3 561 281 87 30 29 13 9 7 2 1 1 0 0 0 0 cohort 1 (n= 409) Recovery location Burt Mullett Crooked Pickerel Black River Recovery year 1 2 3 2 2 0 31 13 9 0 0 0 0 0 0 1 1 0 Burt Lake cohort 2 (n= 687) Recovery year 1 2 3 Recovery location Burt - 62 24 Mullett - 7 1 Crooked - 3 0 Pickerel - 1 1 Black River - 0 0 Mullett Lake cohort 2 (n= 54) Recovery year 1 2 3 Recovery location Burt - 0 0 Mullett - 4 0 Crooked - 0 0 Pickerel - 0 0 Black River - 0 0 79 cohort 3 (n= 2,747) Recovery year 1 2 3 Recovery location Burt - - 122 Mullett - 13 Crooked - 11 Pickerel - 3 Black River - 0 cohort 3 (n= 188) Recovery year 1 2 3 Recovery location Burt - 1 Mullett - 17 Crooked - 0 Pickerel - 0 Black River - 0 Table 1.5: (cont’d) cohort 1 (n= 562) Recovery year 1 2 3 Recovery location Burt 3 2 0 Mullett 0 0 0 Crooked 84 29 11 Pickerel 5 3 0 Black River 0 0 0 cohort 1 (n= 623) Recovery year 1 2 3 Recovery location Burt 5 2 0 Mullett 1 0 0 Crooked 30 4 0 Pickerel 54 26 3 Black River 0 0 0 cohort 1 (n= 261) Recovery year 1 2 3 Recovery location Burt 0 1 0 Mullett 24 8 5 Crooked 0 0 0 Pickerel 0 0 0 Black River 2 0 0 Crooked Lake cohort 2 (n= 529) Recovery year 1 2 3 Recovery location Burt - 4 0 Mullett - 0 0 Crooked - 74 41 Pickerel - 3 3 Black River - 0 0 Pickerel Lake cohort 2 (n= 108) Recovery year 1 2 3 Recovery location Burt - 2 0 Mullett - 0 0 Crooked - 2 2 Pickerel - 10 0 Black River - 0 0 Black River cohort 2 (n= 99) Recovery year 1 2 3 Recovery location Burt - 0 0 Mullett - 6 2 Crooked - 0 0 Pickerel - 0 0 Black River - 0 1 80 cohort 3 (n= 614) Recovery year 1 2 3 Recovery location Burt - 2 Mullett - 0 Crooked - 89 Pickerel - 1 Black River - 0 cohort 3 (n= 326) Recovery year 1 2 3 Recovery location Burt - 3 Mullett - 0 Crooked - 6 Pickerel - 19 Black River - 0 cohort 3 (n= 231) Recovery year 1 2 3 Recovery location Burt - 0 Mullett - 6 Crooked - 0 Pickerel - 0 Black River - 3 APPENDIX 2.0 Supporting tables and figures for Chapter 2: The Influence of Movement Dynamics on Harvest Allocation for Walleye (Sander vitreus) in Intermixed Fisheries in a Chain of Lakes Table 2.1: Location-specific values for parameters used to simulate movement and harvesting dynamics for Michigan’s Inland Waterway. The Abundance (SD) was calculated from a spring mark-recapture study in the Inland Waterway. The estimated fishing mortality, spawning-site fidelity, and post-spawn movement rates and their associated 95% credible intervals were taken from Herbst et al. (Chapter 1). Parameters Abundance Fishing mortality Spawning-site fidelity Post-spawn movement rates Spawning Location Black River Mullett Lake Burt Lake Crooked Lake Pickerel Lake Black River Mullett Lake 477 (54) 2,246 (674) 0.15 (0.02, 0.28) 0.18 (0.13, 0.26) 0.92 (0.85, 0.97) 0.58 (0.38, 0.76) Black River 0.18 (0.04, 0.48) 0.07 (0.00, 0.29) 0.00 (0.00, 0.01) 0.01 (0.00, 0.04) 0.01 (0.00, 0.06) Mullett Lake 0.77 (0.48, 0.92) 0.86 (0.63, 0.94) 0.08 (0.05, 0.10) 0.00 (0.00, 0.02) 0.01 (0.00, 0.04) 81 Burt Lake 19,464 (2,682) 0.19 (0.16, 0.23) 0.99 (0.97, 0.99) Feeding Location Burt Lake 0.03 (0.0, 0.08) 0.06 (0.02, 0.12) 0.90 (0.87, 0.93) 0.06 (0.03, 0.09) 0.09 (0.05, 0.14) Crooked Lake 2,360 (465) 0.31 (0.26, 0.37) 0.96 (0.91, 0.98) Pickerel Lake 4,442 (1,132) 0.28 (0.21, 0.36) 0.80 (0.62, 0.92) Crooked Lake 0.01 (0.0, 0.03) 0.01 (0.0, 0.02) 0.01 (0.0, 0.01) 0.87 (0.81, 0.91) 0.19 (0.13, 0.25) Pickerel Lake 0.01 (0.0, 0.05) 0.01 (0.0, 0.03) 0.01 (0.0, 0.01) 0.06 (0.03, 0.10) 0.70 (0.61, 0.80) Table 2.2: Set of simulation scenarios used to evaluate the implications that movement patterns and differing exploitation rates have on walleye spawning populations within the Inland Waterway. Abundance, movement rates, and recruitment were location-specific and held consistent among all scenarios. Initial abundance was estimated for each location during a spring mark-recapture study in 2011. Natural mortality (M) was held constant at 0.3 and movement rates estimated in Chapter 1 were incorporated for all scenarios. The base scenario (1 and 8) represented perfect implementation of the current harvest control rule (i.e., constant total exploitation rate of 35%) with and without process error, respectively. In scenarios (2-7) movement rates and angling mortality rates (u-angling) were location-specific and taken from estimated posterior distributions (Chapter 1). Scenario u-angling u-spearing Process error 1 0.25 0.1 No 2 Location-specific ua 0.05 No 3 Location-specific ua 0.1 No 4 Location-specific ua 0.2 No 5 Location-specific ua 0.05 Yes 6 Location-specific ua 0.1 Yes 7 Location-specific ua 0.2 Yes 8 0.25 0.1 Yes 82 Table 2.3: Example of the influence that post-spawn movement dynamics have on the total allowable catch (TAC) for each of the locations within the Inland Waterway. Spawning N refers to the estimated abundance during the spring-spawning assessment and the summer N is the abundance at each location after accounting for spearing harvest and post-spawn movements. The spearing TAC was determined using a spearing exploitation rate of 10% and the spawning abundance. The angling TAC was determined using an angling exploitation rate of 25% and the estimate of summer abundance. The angling TAC adjustments were determined by evaluating the difference between the two movement scenarios, with negative values indicating that fewer fish would need to be harvested to account for movement dynamics while maintaining the u = 0.35 for the spawning populations throughout the waterway. Perfect Implementation of Harvest Control Rule (u = 0.35) Movement Ignored Location Burt Mullett Crooked Pickerel Black River Totals Spawning Spearing Summer Angling N TAC N TAC 19,464 1,946 17,518 4,380 2,246 225 2,021 505 2,360 236 2,124 531 4,442 444 3,998 1000 477 48 429 107 28,989 2,899 26,090 6,523 Angling TAC adjustment accounting for movement Spawning Spearing Summer Angling so u = 0.35 for each N TAC N TAC spawning population 19,464 1,946 16,387 4,097 -283 2,246 225 3,510 878 372 2,360 236 2,807 702 171 4,442 444 3,105 776 -224 477 48 281 71 -35 28,989 2,899 26,090 6,523 0 Movement Incorporated 83 Table 2.4: Mean exploitation rates for each spawning population under different levels of uncertainty and spearing exploitation (us= 0.05, 0.10, and 0.20) and observed lake-specific angling mortality rates (ua). The uncertainty evaluated was process error on total instantaneous mortality (Z) over summer within the population dynamics model. (No error = process error excluded; Error = process error included; Diff = difference between the two approaches). Low (us = 0.05) Medium (us = 0.10) High (us = 0.20) Location No error Error Diff No error Error Diff No error Error Diff Burt Lake 0.20 0.20 0.00 0.24 0.24 0.00 0.32 0.32 0.00 Mullett Lake 0.18 0.19 -0.01 0.23 0.24 -0.01 0.31 0.33 -0.02 Crooked Lake 0.27 0.26 0.01 0.31 0.30 0.01 0.38 0.39 -0.01 Pickerel Lake 0.25 0.25 0.00 0.29 0.29 0.00 0.37 0.37 0.00 Black River 0.18 0.18 0.00 0.22 0.22 0.00 0.26 0.31 -0.05 84 Figure 2.1: Map of northern Michigan’s Inland Waterway that consists of four lakes (Burt, Crooked, Mullett, and Pickerel) and four major connecting rivers (north to south through the lakes: Cheboygan River, Black River, Indian River, and Crooked River). 85 Fish exhibiting spawning site fidelity Figure 2.2: Conceptual diagram depicting the modeling process of how a spawning population (e.g., Burt Lake population) is tracked and projected through time in the Inland Waterway using the stochastic simulation model. Abundance is initially estimated using mark-recapture methods during spawning. Following the assessment of abundance ( ̂) the population is subjected to tribal harvest (i.e., spearing (us)) within the spawning grounds. After spawning, lake-specific spawning populations exhibit postspawn movement (ϕ) and are subjected to total instantaneous mortality (Z = angling (ua) + natural (M)) in summer feeding locations. The fraction of the spawning population that survives during time t then returns to spawning grounds (i.e., spawning-site fidelity (ψ)) or remains in the location that they resided in during summer feeding. New additions represent immigrants from other spawning populations that do not exhibit ψ. During time t+1 the spawning populations are projected forward with the addition of constant recruitment that is specified at a level consistent with spawning population-specific harvest from time t. Locations abbreviations: BL = Burt Lake, ML = Mullett Lake, CL = Crooked Lake, PL = Pickerel Lake, and BR = Black River. 86 Figure 2.3: Exploitation rates for each spawning population under the base scenarios where the harvest control rule was implemented without error (panel A; u = 0.35). Panel B represents the base scenario, but also includes a process error term on total instantaneous mortality (Z) within the population dynamics model. The red line indicates the target exploitation rate set forth as the waterway’s harvest control rule (u = 0.35). 87 Figure 2.4: Influence of varying levels of spearing exploitation rates (Scenario 2: us= 0.05 (low); Scenario 3: us = 0.10 (medium); Scenario 4: us = 0.20 (high)) and location-specific observed fishing mortality on exploitation rate (u) for each spawning population within the waterway. The exploitation rates under the base scenario that represents perfect implementation of the harvest control rule (HCR) is also provided for reference. The red line indicates the 95% CI for the maximum exploitation rate (u = 0.35), given the spawning population abundance estimate and standard deviation. 88 APPENDIX 3.0 Supporting tables and figures for Chapter 3: Walleye Foraging Ecology in an Interconnected Chain of Lakes Influenced by NonNative Species. Table 3.1: Characteristics of lakes in the Inland Waterway. The percent littoral is defined as lake area ≤ 4.6 m total depth. All values are based off the 2013 mean. Total Phosphorus (TP) represents the bottom-to-surface mean. Crooked Lake Secchi and Chl-a are from 2011, which were the most recent data. Water chemistry data were provided by Tip of the Mitt Watershed Council (Kevin Cronk personal communication). Lake size and bathymetric information are unpublished MDNR data. Lake name Surface area (km2) Percent littoral* Mean depth (m) Max depth (m) Secchi TP Chl-a (m) (ug/L) (ug/L) Non-natives present Burt 70.4 28 7.8 22.3 5.7 3.5 1.23 Round Gobies, zebra mussels, alewife Crooked 9.9 80.3 2.7 20.7 2.9 12.7 3.04 Zebra mussels Mullett 67.3 31.1 10.4 45.1 5.2 2.5 0.94 Round Gobies, zebra mussels, alewife, rainbow smelt Pickerel 4.4 78.8 3.2 20.4 3.8 2.8 1.54 Zebra mussels 89 Table 3.2: Catch-per-unit-effort (CPUE: fish/net) of fish in forage gill nets fished during May-August, 2011 and June-August, 2012 in the Inland Waterway. Only species that represent potential prey items for walleye and other predatory fish are listed, so incidental catches of predator species were omitted. Species Yellow perch Spottail shiner Alewife Trout perch Cisco Rainbow smelt Round goby Totals Mullett Lake 2011 2012 12.6 8.3 0.4 0.5 0.5 2.6 0.1 0.2 0.0 0.0 0.2 0.3 0.1 0.0 13.9 11.9 Burt Lake 2011 2012 6.7 12.6 0.5 0.9 0.1 0.1 0.1 0.0 0.3 0.2 0.0 0.0 0.0 0.0 7.7 13.8 90 Crooked Lake 2011 2012 5.6 9.3 0.8 0.3 0.1 0.1 0.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.6 9.7 Pickerel Lake 2011 2012 2.0 2.4 0.6 0.6 0.0 0.0 0.2 0.1 0.0 0.0 0.0 0.0 0.0 0.0 2.8 3.1 Table 3.3: Seasonal diet composition (% wet weight (g)) for walleye from the Inland Waterway, 2011-2013. Seasons were defined as spring (April-June), summer (July-September), and fall/winter (October-March). N is the total number of stomachs that were analyzed for prey items. Unk. fish category consisted of unidentifiable bones and the “other” category consisted of miscellaneous fish species such as emerald and spottail shiners. Diet item Location Aquatic insects Yellow perch Round goby Burt Lake Crooked Lake Mullett Lake Pickerel Lake 62.5 58.7 24.6 53.3 10.3 2.8 62.3 0.0 6.6 0.0 0.4 0.0 Burt Lake Crooked Lake Mullett Lake Pickerel Lake 10.3 12.6 3.8 15.1 21.1 13.3 38.8 9.8 36.4 0.3 17.2 0.0 Burt Lake Crooked Lake Mullett Lake 0.3 0.8 0.0 33.3 7.0 10.2 42.9 0.0 77.3 Pickerel Lake - - - Crayfish Unk. fish Spring 17.3 2.0 6.3 29.1 0.0 4.7 43.6 1.3 Summer 22.1 7.3 52.1 7.5 6.9 19.3 54.8 12.8 Fall/winter 5.9 15.4 0.0 84.7 0.0 0.1 - 91 - Alewife and smelt Other N % empty 0.0 0.0 0.0 0.0 1.3 3.1 8.0 1.8 579 46 32 60 25.4 21.7 40.6 38.3 0.0 0.0 6.9 0.0 2.8 14.2 7.1 7.5 148 191 143 67 50.7 46.1 55.2 47.8 0.0 0.0 2.9 2.2 7.5 9.5 185 23 10 27.0 30.4 10.0 - - 0 - Table 3.4: Stable isotope sample sizes (N), estimates of trophic position and percent littoral (SD) along with the standard ellipse area corrected for sample size (SEAc) for walleye within the Inland Waterway, in northern Michigan. Trophic position Percent littoral lower 95% SEAc Mean upper 95% Lake N Mullett 20 4.36 (0.29) 86.4 (12.63) 6.4 7.3 7.7 Burt 19 4.55 (0.14) 72.8 (6.24) 2.5 2.8 2.9 Crooked 19 4.55 (0.13) 96.2 (5.42) 1.5 1.8 1.9 Pickerel 20 4.39 (0.10) 68.1 (16.91) 4.1 4.8 5.1 92 Figure 3.1: Map of northern Michigan’s Inland Waterway that consists of four lakes (Burt, Crooked, Mullett, and Pickerel) and four major connecting rivers (north to south through the lakes: Cheboygan River, Black River, Indian River, and Crooked River). 93 Figure 3.2: Stable isotope bi-plots illustrating mean trophic position and percent littoral (± 1 SE) for multiple species within the communities of Mullett, Burt, Crooked, and Pickerel lakes within the Inland Waterway. Species-specific abbreviations include: ZOP = zooplankton, ALE = alewife, YEP = yellow perch, WAE = walleye, RGB = round goby, SPO = spottail shiner, CRAY = crayfish species, GST = gastropods, HEX= Hexagenia spp., and SMT = rainbow smelt. 94 Figure 3.3: Stable isotope bi-plot representing the percent littoral and trophic position for walleye collected from each of the lakes within the Inland Waterway, MI. The lines enclose the standard ellipse areas (SEAc) and represent the total isotopic niche area for walleye within each of the lakes. The symbols are data points that indicate the isotopic signature of individual walleye from Burt (open triangle), Crooked (cross hatch), Mullett (open circle), and Pickerel (“x”) lakes. 95 APPENDIX 3.1 Stable Isotope summary statistics (δ13C, δ15N, trophic position (TP), and percent littoral). Table 3.5: Stable Isotope summary statistics (δ13C, δ15N, trophic position (TP), and percent littoral) for each species analyzed by lake. N indicates the total sample size collected, average length was measured in TL (nearest mm), values in parentheses represent the SD, covar is the covariance between carbon and nitrogen values, and * indicates composite samples. Mullett Species Alewife Crayfish Gastropod Hexagenia Rainbow smelt Round goby Spottail shiner Walleye Yellow Perch Zooplankton N Length 20 157 (20.9) 20 28 (4.4) 5* 0 12 147 (15.8) 20 64 (20.7) 0 20 450 (110.3) 20 100 (14.4) 5* - Burt δ13C δ15N Covar -32.0 (0.64) 9.1 (0.66) -0.213 -26.9 (0.37) 7.4 (0.54) 0.081 -27.3 (0.21) 5.9 (0.14) 0.012 -29.8 (0.24) 9.9 (0.43) 0.039 -27.3 (0.72) 9.8 (0.44) 0.076 -27.6 (0.51) 12.5 (0.98) 0.381 -28.5 (1.04) 8.7 (0.47) 0.053 -30.6 (1.54) 5.7 (0.83) 0.954 96 N 1 20 5* 3* 0 20 20 20 20 5 Length 171 (-) 35 (3.0) 66 (16.7) 74 (18.3) 408 (43.3) 108 (26.1) - δ13C -32.0 (-) -26.6 (1.08) -26.4 (0.33) -31.7 (1.14) -28.1 (1.18) -27.6 (1.12) -28.3 (0.43) -30.3 (1.68) -33.4 (0.25) δ15N 10.2 (-) 8.3 (0.61) 7.6 (0.15) 8.6 (0.79) 10.8 (0.50) 10.2 (1.18) 13.7 (0.47) 10.7 (1.01) 7.3 (0.15) Covar -0.529 0.001 -0.020 -0.078 -0.889 0.012 1.403 0.025 Table 3.5: (cont’d) Species Alewife Crayfish Gastropod Hexagenia Rainbow smelt Round goby Spottail shiner Walleye Yellow Perch Zooplankton N 3 20 5* 4 0 0 20 20 20 5* Length 86 (2.0) 34 (6.2) 92 (9.3) 411 (44.5) 99 (13.3) - Crooked δ13C -29.9 (0.42) -27.5 (0.84) -28.3 (0.69) -32.8 (0.51) -30.9 (0.73) -28.3 (0.55) -30.3 (0.71) -34.1 (0.54) 15 δ N 10.9 (0.10) 8.8 (0.83) 8.4 (0.28) 7.1 (0.51) 11.3 (0.54) 14.3 (0.43) 11.1 (0.56) 8.9 (0.20) Covar -0.027 -0.530 -0.020 0.159 0.004 0.102 -0.098 -0.012 97 N 0 17 5* 0 0 0 20 20 20 5* Length 33 (6.4) 87 (12.5) 380 (47.2) 94 (10.0) - Pickerel δ13C δ15N -29.9 (0.87) 8.9 (1.17) -28.4 (1.03) 7.4 (0.21) -32.1 (1.29) 10.1 (0.59) -30.0 (0.97) 12.9 (1.04) -31.4 (1.35) 10.5 (0.50) -33.6 (0.56) 8.2 (0.15) Covar 0.025 0.066 -0.129 0.076 0.314 0.012 Table 3.5: (cont’d) Species Alewife N 20 Length 157 (20.9) Mullett % Littoral 0.5 (2.2) TP Covar 2.98 (0.19) -0.194 N 1 Length 171 (-) Burt % Littoral 19.6 (-) TP 2.8 (-) Covar - Crayfish Gastropod 20 5* 28 (4.4) - 99.1 (2.4) 97.5 (3.1) 2.5 (0.16) 2.1 (0.05) 0.021 0.083 18 5* 35 (2.8) - 93.8 (11.9) 97.5 (5.7) 2.3 (0.17) 2.7 (0.05) -1.538 0.127 Hexagenia 0 - - - - 3* - 0.73 (1.3) 3.1 (0.25) 0.256 Rainbow smelt 11 147 (15.8) 23.4 (7.2) 3.2 (0.13) 0.339 0 - - - - Round goby 20 64 (20.7) 89.8 (10.5) 3.2 (0.13) 0.258 20 66 (16.7) 75.2 (17.0) 3.0 (0.16) -0.131 Spottail shiner Walleye Yellow Perch 0 20 20 450 (110.3) 100 (14.4) 86.4 (12.6) 59.8 (27.4) 4.0 (0.29) 2.9 (0.13) 2.915 0.299 20 19 19 74 (18.3) 410 (43.0) 106 (25.7) 82.8 (15.9) 72.8 (6.2) 44.3 (22.3) 2.9 (0.35) 3.9 (0.15) 2.9 (0.30) -3.832 0.109 5.726 Zooplankton 4* - 16.5 (33.1) 2.0 (0.27) 8.813 5 - 0.0 (0.0) 2.7 (0.05) 0.000 98 Table 3.5: (cont’d) Species Alewife Crayfish Gastropod Hexagenia Rainbow smelt Round goby Spottail shiner Walleye Yellow Perch Zooplankton N 3 20 4* 4 0 0 17 19 19 5* Length 86 (2.0) 34 (6.2) 92 (9.7) 412 (45.6) 99 (13.4) - Crooked % Littoral 71.6 (7.1) 98.2 (5.7) 84.8 (18.8) 0.0 (0.0) 53.6 (12.5) 96.2 (5.4) 65.2 (12.1) 0.0 (0.0) TP 3.5 (0.06) 2.9 (0.24) 2.8 (0.10) 2.5 (0.17) 3.7 (0.15) 4.6 (0.13) 3.6 (0.17) 3.0 (0.08) Covar -0.400 -1.032 -0.242 0.000 0.107 0.375 -0.549 0.000 99 N 0 17 5* 0 0 0 20 20 20 5* Length 33 (6.4) 87 (12.5) 380 (47.2) 94 (10.0) - Pickerel % Littoral 69.9 (16.8) 91.5 (12.3) 30.0 (24.0) 68.1 (16.9) 44.8 (21.7) 4.2 (7.5) TP 2.2 (0.36) 1.8 (0.05) 2.6 (0.16) 3.5 (0.09) 2.7 (0.14) 2.0 (0.04) Covar 0.238 -0.038 -0.769 0.391 1.276 -0.011 APPENDIX 3.2 Analysis of covariance results to determine the influence of total length on trophic position for walleye within the Inland Waterway. Figure 3.4: Plot of analysis of covariance results for the influence of total length on trophic position for each lake. 100 Table 3.6: Trophic position by length regression line specifics along with the overall mean length adjustment to determine biological influence associated with the significant interaction term that indicated trophic position was influenced by total length differently by location Lake Mullett Burt Crooked Pickerel N 20 19 19 20 Length (mm) 450 (110.3) 410 (43.0) 412 (45.6) 380 (47.2) Slope 0.0023 0.0002 0.0018 0.0009 y-intercept 2.97 3.79 3.81 3.15 Observed mean trophic position 4.0 3.9 4.6 3.5 101 Adjusted trophic position for length = 412 mm 3.9 3.9 4.6 3.5 REFERENCES 102 REFERENCES Aires-da-Silva, A.M., M.N. Maunder, V.F. Gallucci, N.E. Kohler, and J.J. Hoey. 2009. A spatially structured tagging model to estimate movement and fishing mortality rates for blue shark (Prionace glauca) in the North Atlantic Ocean. Marine and Freshwater Research 60:1029-1043. Ames, E.P. 2004. Atlantic cod stock structure in the Gulf of Maine. Fisheries 29:10-28. Anganuzzi, A., R. Hilborn, and J.R. Skalski. 1994. Estimation of size selectivity and movement rates from mark-recovery data. Canadian Journal of Fisheries and Aquatic Sciences 51:734-742. Beisner, B.E., A.R. Ives, and S.R. Carpenter. 2003. The effects of an exotic fish invasion on the prey communities of two lakes. Journal of Animal Ecology 72:331-342. Benson, A. J. 2014. Zebra mussel sightings distribution. Retrieved 10/13/2014 from http://nas.er.usgs.gov/taxgroup/mollusks/zebramussel/zebramusseldistribution.aspx. Berger, A.M., M.L. Jones, Y. Zhao, and J.R. Bence. 2012. Accounting for spatial population structure at scales relevant to life history improves stock assessment: the case for Lake Erie walleye Sander vitreus. Fisheries Research 115-116:44-59. Beverton, R. J. H., and S.J. Holt.1957. On the dynamics of exploited fish populations, Fishery Investigations Series II Volume XIX, Ministry of Agriculture, Fisheries and Food. Billington, N., C.C. Wilson, and B.L. Sloss. 2011. Distribution and population genetics of walleye and sauger. Pages 105-132 in B.A. Barton, editor. Biology, management, and culture of walleye and sauger. American Fisheries Society, Bethesda, Maryland. Bjorkvoll, E., V. Grotan, S. Aanes, B.E. Saether, S. Engen, and R. Aanes. 2012. Population dynamics and life-history variation in marine fish species. The American Naturalist 180:372-387. Booth, A.J. 2000. Incorporating the spatial component of fisheries data into stock assessment models. ICES Journal of Marine Science 57:858-865. Brooks, J.L. and S.I. Dodson. 1965. Predation, body size, and composition of plankton. Science 150:28-35. Brownie, C., J.E. Hines, J.D. Nichols, K.H. Pollock, and J.B. Hestbeck. 1993. Capture-recapture studies for multiple strata including non-Markovian transitions. Biometrics 49:11731187. 103 Buckland, S.T., I.B.J. Goudie, and D.L. Borchers. 2000. Wildlife population assessment: past developments and future directions. Biometrics 56:1-12. Buckland, S.T., K.B. Newman, C. Fernandez, L. Thomas, and J. Harwood. 2007. Embedding population dynamics models in inference. Statistical Science 22:44-58. Butterworth, D.S., and A.E. Punt. 1999. Experiences in the evaluation and implementation of management procedures. ICES Journal of Marine Science 56:985-998. Campbell, L.M., R. Thacker, D. Barton, D.C.G. Muir, D. Greenwood, and R.E. Hecky. 2009. Re-engineering the eastern Lake Erie littoral food web: the trophic function of nonindigenous Ponto-Caspian species. Journal of Great Lakes Research 35:224-231. Chipps, S. R. and B. D. S. Graeb. 2011. Feeding ecology and energetics. Pages 303-319 in B. A. Barton, editor. Biology, management, and culture of walleye and sauger. American Fisheries Society, Bethesda, Maryland. Clark, J. 2007. Models for ecological data: an introduction. Princeton University Press. 152pp. Conroy, M.J., and W.W. Blandin. 1984. Geographic and temporal differences in band reporting rates for American black ducks. Journal of Wildlife Management 48:23-36. Cope, J.M., and A.E. Punt. 2009. Drawing the lines: resolving fishery management unit with simple fisheries data. Canadian Journal of Fisheries and Aquatic Sciences 69:1256-1273. Crowe, W. R. 1962. Homing behavior in walleyes. Transactions of the American Fisheries Society 91:350-354. Cucherousset, J. and J.D. Olden. 2011. Ecological impacts of nonnative freshwater fishes. Fisheries 36:215-230. DePhilip, M. M., J.S. Diana, and D. Smith. 2005. Movement of walleye in an impounded reach of the Au Sable River, Michigan, USA. Environmental Biology of Fishes 72:455-463. Deroba, J.J., and J.R. Bence. 2008. A review of harvest policies: understanding relative performance of control rules. Fisheries Research 94:210-223. Dietrich, J.P., B.J. Morrison, and J.A. Hoyle. 2006. Alternative ecological pathways in the eastern Lake Ontario food web: Round Goby in the diet of Lake Trout. Journal of Great Lakes Research 32:395-400. Dupuis, J.A. 1995. Bayesian estimation of movement and survival probabilities from capturerecapture data. Biometrika 82: 761-772. Dupuis, J.A., and C.J. Schwarz. 2007. A Bayesian approach to the multistate Jolly-Seber capture-recapture model. Biometrics 63: 1015-1022. 104 Ellison, A.M. 2004. Bayesian inference in ecology. Ecology Letters 7:509-520. Fielder, D. G., and M. V. Thomas. 2006. Fish population dynamics of Saginaw Bay, Lake Huron 1998 – 2004. Michigan Department of Natural Resources, Fisheries Research Report. No. 2083. Ann Arbor. Forney, J. L. 1974. Interactions between yellow perch abundance, walleye predation, and survival of alternate prey in Oneida Lake, New York. Transactions of the American Fisheries Society 103:15-24. France, R. 1995. Differentiation between littoral and pelagic food webs in lakes using stable carbon isotopes. Limnology and Oceanography 40:1310-1313. Frederick, S.W., and R.M. Peterman. 1995. Choosing fisheries harvest policies: when does uncertainty matter? Canadian Journal of Fisheries and Aquatic Sciences 52:291-306. Fu, C., and L.P. Fanning. 2004. Spatial considerations in the management of Atlantic cod off Nova Scotia, Canada. North American Journal of Fisheries Management 24:775-784. Fuiman, L.A., and A.E. Magurran. 1994. Development of predator defenses in fishes. Reviews of Fish Biology and Fisheries 4:145-183. Garcia-Berthou, E. 2007. The characteristics of invasive fishes: what has been learned so far? Journal of Fish Biology 71:33-55. Gelman, A., and D.R. Rubin. 1992. A single series from the Gibbs sampler provides a false sense of security. Pages 625-631 in J.M. Bernardo, J.O. Berger, A.P. Dawid, and A.F.M. Smith, editors. Bayesian statistics 4. Oxford University Press, Oxford, U.K. Gelman, A., J. Carlin, H.S. Stern, and D.B. Rubin. 2004. Bayesian data analysis. Chapman and Hall/CRC Press, Boca Raton, Florida. Gimenez, O., V. Rossi, R. Choquet, C. Dehais, B. Doris, H. Varella, J. Vila, and R. Pradel. 2007. State-space modelling of data on marked individuals. Ecological Modelling 206:431-438. Goethel, D.R., T.J. Quinn, and S.X. Cadrin. 2011. Incorporating spatial structure in stock assessment: movement modeling in marine fish population dynamics. Reviews in Fisheries Science 19:119-136. Guan, W., J. Cao, Y. Chen, and M. Cieri. 2013. Impacts of population and fishery spatial structures on fishery stock assessment. Canadian Journal of Fisheries and Aquatic Sciences 70:1178-1189. Guzzo, M.M., G.D. Haffner, N.D. Legler, S.A. Rush, and A.T. Fisk. 2013. Fifty years later: trophic ecology and niche overlap of a native and non-indigenous fish species in the western basin of Lake Erie. Biological Invasions 15:1695-1711. 105 Hanchin, P. A., R. D. Clark, R. N. Lockwood, and N. A. Godby. 2005a. The fish community and fishery of Crooked and Pickerel lakes, Emmet County, Michigan with emphasis on walleyes and northern pike. Michigan Department of Natural Resources, Fisheries Special Report 34, Ann Arbor. Hanchin, P. A., R. D. Clark, R. N. Lockwood, and T. A. Cwalinski. 2005b. The fish community and fishery of Burt Lake, Cheboygan County, Michigan in 2001-02 with emphasis on walleyes and northern pike. Michigan Department of Natural Resources, Fisheries Special Report 36, Ann Arbor. Hansen, M.J. 1989. A walleye population model for setting harvest quotas. Wisconsin Department of Natural Resources, Bureau of Fisheries Management, Fish Management Report 143, Madison, WI. Hansen, M.J., A.H. Fayram, and S.P. Newman. 2011. Natural mortality in relation to age and fishing mortality on walleyes in Escanaba Lake, Wisconsin, during 1956-2009. North American Journal of Fisheries Management 31:506-514. Hayes, D. B., J. R. Bence, T. J. Kwak, and B. E. Thompson. 2007. Abundance, biomass, and production. Pages 327-374 in C. S. Guy and M. L. Brown, editors. Analysis and interpretation of freshwater fisheries data. American Fisheries Society, Bethesda, Maryland. Hayes, D.B., C.P. Ferreri, and W.W. Taylor. 2012. Active Fish Capture Techniques. Pages 267304 in A.V. Zale, D.L. Parrish, and T.M. Sutton, editors. Fisheries Techniques, 3rd edition. American Fisheries Society, Bethesda, MD. Hayes, D, B. Burroughs, and B. Thompson. 2014. Resource Management in the Face of Uncertainty. Pages 267-271 in W. Taylor, A. Lynch, N. Leonard, editors. Future of Fisheries: Perspectives for Emerging Professionals. American Fisheries Society, Bethesda, MD. He, J.X., J.R. Bence, C.P. Madenjian, S.A. Pothoven, N.E. Dobiesz, D.G. Fielder, J.E. Johnson, M.P. Ebener, R.A. Cottrill, L.C. Mohr, and S.R. Koproski. 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 Science 72:7-23. Hendrix, A.N., J. Straley, C.M. Gabriele, and S.M. Gende. 2012. Bayesian estimation of humpback whale (Megaptera novaeangliae) population abundance and movement patterns in southeastern Alaska. Canadian Journal of Fisheries and Aquatic Science 69:1783-1797. Henny C.J., and K.P. Burnham. 1976. A reward band study of mallards to estimate band reporting rates. Journal of Wildlife Management 40:1-14. 106 Herbst, S.J., J.E. Marsden, and B.F. Lantry. 2013. Lake whitefish diet, condition, and energy density in Lake Champlain and the lower four Great Lakes following dreissenid invasions. Transactions of the American Fisheries Society 142:388-398. Hilborn, R. 1990. Determination of fish movement patterns from tag recoveries using maximum likelihood estimators. Canadian Journal of Fisheries and Aquatic Science 47:635-643. Hobson, K.A., A. Ofukany, D.X. Soto, L.I. Wassenaar. 2012. An isotopic baseline (δ13C, δ15N) for fishes of Lake Winnipeg: implications for investigating impacts of eutrophication and invasive species. Journal of Great Lakes Research 38:58-65. Holeck, K.T., E.L. Mills, H.J. MacIsaac, M.R. Dochoda, R.I. Colautti, and A. Ricciardi. 2004. Bridging troubled waters: understanding links between biological invasions, transoceanic shipping, and other entry vectors in the Laurentian Great Lakes. Bioscience 54:919-929. Hutchinson, W.F. 2008. The dangers of ignoring stock complexity in fishery management: the case of the North Sea cod. Biology Letters 4:693-695. Isermann, D.A. and C.T. Knight. 2005. Potential effects of jaw tag loss on exploitation estimates for Lake Erie walleyes. North American Journal of Fisheries Management 25:557-562. Jackson, A.L., R. Inger, A.C. Parnell, and S. Bearhop. 2011. Comparing isotopic niche widths among and within communities: SIBER – Stable Isotope Bayesian Ellipses in R. Journal of Animal Ecology 80:595-602. Jackson, M.C., I. Donohue, A.L. Jackson, J.R. Britton, D.M. Harper, and J. Grey. 2012. Population-level metrics of trophic structure based on stable isotopes and their application to invasion ecology. Plos ONE 7:e31757. Johnson, T.B., D.B. Bunnell, and C.T. Knight. 2005. A potential new energy pathway in central Lake Erie: the round goby connection. Journal of Great Lakes Research 31:238-251. Jones, M.S., J.P. Goettl, and S.A. Flickinger. 1994. Changes in walleye food habits and growth following a rainbow smelt introduction. North Amercian Journal of Fisheries Management 14:409-414. Kaemingk, M.A., T.L. Galarowicz, J.A. Clevenger, D.A. Clapp, and H.L. Lenon. 2012. Fish assemblage shifts and population dynamics of smallmouth bass in the Beaver Archipelago, northern Lake Michigan: a comparison between historical and recent time periods amidst ecosystem changes. Transactions of the American Fisheries Society 141:550-559. Katsukawa, T. 2004. Numerical investigation of the optimal control rule for decision-making in fisheries management. Fisheries Science 70:123-131. 107 Kéry, M., and M. Schaub. 2011. Bayesian Population Analysis Using WinBUGS - a Hierarchical Perspective. Academic Press. Kilian, J.V., R.J. Klauda, S. Widman, M. Kashiwagi, R. Bourquin, S. Weglein, and J. Schuster. 2012. An assessment of a bait industry and angler behavior as a vector of invasive species. Biological Invasions 14:1469-1481. King, R. 2014. Statistical ecology. Annual Review of Statistics and its Applications 1:401-426. Koenigs, R.P., R.M. Bruch, and K.K. Kamke. 2013. Impacts of anchor tag loss on walleye management in the Winnebago system. North American Journal of Fisheries Management 33:909-916. Kornis, M.S., N. Mercado-Silva, and M.J. Vander Zanden. 2012. Twenty years of invasion: a review of round goby Neogobius melanostomous biology, spread and ecological implications. Journal of Fish Biology 80:235-285. Kornis, M.S., S. Sharma, and M.J. Vander Zanden. 2013. Invasion success and impact of an invasive fish, round goby, in Great Lakes tributaries. Diversity and Distributions 19:184198. Krueger, D.M. and T.R. Hrabik. 2005. Food web alteration that promote native species: the recovery of cisco (Coregonus artedi) populations through management of native piscivores. Canadian Journal of Fisheries and Aquatic Science 62: 2177-2188. Layman, C.A., D.A. Arrington, C.G. Montana, and D.M. Post. 2007. Can stable isotope ratios provide for community-wide measures of trophic structure? Ecology 88:42-48. Li, Y., J.R. Bence, and T.O. Brenden. 2014. An evaluation of alternative assessment approaches for intermixing fish populations: a case study with Great Lakes lake whitefish. ICES Journal of Marine Science. In press. Litvak, M.K., and N.E. Mandrak. 1993. Ecology of freshwater bait-fish use in Canada and the United States. Fisheries 18:6-13. Lodge, D.M. 1993. Biological invasions: lessons for ecology. Trends in Ecology and Evolution 8:133-137. Lumb, C.E., T.B. Johnson, H.A. Cook, and J.A. Hoyle. 2007. Comparison of Lake Whitefish (Corgeonus clupeaformis) growth, condition, and energy density between Lakes Erie and Ontario. Journal of Great Lakes Research 33:314-325. Madenjian, C.P., M.A. Stapanian, L.D. Witzel, D.W. Einhouse, S.A. Pothoven, and H.L. Whitford. 2011. Evidence for predatory control of the invasive round goby. Biological Invasions 13:987-1002. 108 Madenjian, C.P., R. O’Gorman, D.B. Bunnell, R.L. Argyle, E.F. Roseman, D.M. Warner, J.D. Stockwell, and M.A. Stapanian. 2013. Adverse effects of alewives on Laurentian Great Lakes fish communities. North American Journal of Fisheries Management 28:263-282. Meng, X.I. 1994. Posterior predictive p-values. Ann. Stat. 22:1142-1160. Mercado-Silva, N., G.G. Sass, B.M. Roth, S. Gilbert, and M.J. Vander Zanden. 2007. Impact of rainbow smelt (Osmerus mordax) invasion on walleye (Sander vitreus) recruitment in Wisconsin lakes. Canadian Journal of Fisheries and Aquatic Sciences 64:1543-1550. Minagawa, M. and E. Wada. 1984. Stepwise enrichment of 15N along food chains: Further evidence and the relation between δ15N and animal age. Geochimica et Cosmochimica Acta 48:1135-1140. Molton, K.J., T,O. Brenden, and J.R. Bence. 2012. Control rule performance for intermixing lake whitefish populations in the 1836 Treaty waters of the Great Lakes: a simulation-based evaluation. Journal of Great Lakes Research 38:686-698. Molton, K.J., T.O. Brenden, and J.B. Bence. 2013. Harvest levels that conserve spawning biomass can provide larger and more stable and sustainable yields in intermixed fisheries. Fisheries Research 147:264-283. Morishima, G.S., and K.A. Henry. 1999. The history and status of Pacific Northwest Chinook and coho salmon ocean fisheries and prospects for sustainability. In Sustainable Fisheries Management: Pacific Salmon, pp. 219-236. Ed. By E.E. Knudsen, C.R. Steward, D.R. MacDonald, J.E. Williams, and D.W. Reiser. CRC Press, Boca Raton, FL. 752pp. Moyle, P.B., and J.J. Cech. 2004. Fishes – An Introduction to Ichthyology fifth edition. Pearson Benjamin Cummings, San Francisco, C.A. Nate, N.A., M.A. Bozek, M.J. Hansen, and S.W. Hewett. 2000. Variation in walleye abundance with lake size and recruitment source. North American Journal of Fisheries Management 20:119-126. Nieland, J.L., M.J. Hansen, M.J. Seider, and J.J. Deroba. 2008. Modeling the sustainability of lake trout fisheries in eastern Wisconsin waters of Lake Superior. Fisheries Research 94:304-314. Olson, D. E., and W.J. Scidmore. 1962. Homing behavior of spawning walleye. Transactions of the American Fisheries Society 91:355-361. Overman, N. C. and D. L. Parrish. 2001. Stable isotope composition of walleye: 15N accumulation with age and area-specific differences in δ13C. Canadian Journal of Fisheries and Aquatic Sciences 58:1253-1260. 109 Perga, M.E. and D. Gerdeaux. 2004. Changes in the δ13C of pelagic food webs: the influence of lake area and trophic status on the isotopic signature of whitefish (Coregonus lavaretus). Canadian Journal of Fisheries and Aquatic Sciences 61:1485-1492. Pine, W.E., K.H. Pollock, J.E. Hightower, T.J. Kwak, and J.A. Rice. 2003. A review of tagging methods for estimating fish population size and components of mortality. Fisheries Research 28:10-23. Pollock, K.H., J.M. Hoenig, and C.M. Jones. 1991. Estimation of fishing and natural mortality when a tagging study is combined with creel or port sampling. Pages 423-434 in D. Guthrie, J.M. Hoenig, M. Holliday, C.M. Jones, M.J. Mills, S.A. Moberly, K.H. Pollock, and D.R. Talhelm, editors. Creel and angler surveys in fisheries management. American Fisheries Society, Symposium 12, Bethesda, Maryland. Porch, C.E., S.C. Turner, and J.E. Powers. 2001. Virtual population analyses of Atlantic Bluefin tuna with alternative models of transatlantic migration: 1970-1997. ICCAT Collected Volume of Scientific Papers 52:1022-1045. Post, D.M. 2002. Using stable isotopes to estimate trophic position: models, methods, and assumptions. Ecology 83:703-718. Pothoven, S.A., and C.P. Madenjian. 2008. Changes in consumption by Alewives and Lake Whitefish after dreissenid mussel invasions in Lakes Michigan and Huron. North American Journal of Fisheries Management 28:308-320. Pothoven, S.A., and C.P. Madenjian. 2013. Increased piscivory by lake whitefish in Lake Huron. North American Journal of Fisheries Management 33:1194-1202. Pothoven, S.A., and T.F. Nalepa. 2006. Feeding ecology of Lake Whitefish in Lake Huron. Journal of Great Lakes Research 32:489-501. Punt, A.E. 2006. The FAO precautionary approach after almost 10 years: have we progressed towards implementing simulation-tested feedback-control management systems for fisheries management? Natural Resource Modeling 19:441-464. Quinn, T.J., and R.B. Deriso. 1999. Quantitative Fish Dyanmics. Oxford UniversityPress, New York, N.Y. R Core Team. 2010. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0. Available from http://www.R-project.org/. R Development Core Team. 2012. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna. 110 Rasmussen, P. W., D.M. Heisey, S.J. Gilbert, R.M. King, and S.W. Hewett. 2002. Estimating postspawning movement of walleyes among interconnected lakes of northern Wisconsin. Transactions of the American Fisheries Society 131:1020-1032. Reeves, K. A. In review. Mille Lacs forage study, 2006-2008. Minnesota Department of Natural Resources Section of Fisheries, Investigational Report. Reyjol, Y., P. Brodeur, Y. Mailhot, M. Mingelbier, and P. Dumont. 2010. Do native predators feed on non-native prey? The case of round goby in a fluvial piscivorous fish assemblage. Journal of Great Lakes Research 36:618–624. Ricciardi, A., M.F. Hoopes, M.P. Marchetti, and J.L. Lockwood. 2013. Progress toward understanding the ecological impacts of nonnative species. Ecological Monographs 83:263-282. Rice, J. A. 2007. Mathematical statistics and data analysis. 3rd edition. Thompson Higher Education, Belmont, California, USA. Roell, M.J., and D.J. Orth. 1993. Trophic basis of production of stream-dwelling smallmouth bass, rock bass, and flathead catfish in relation to invertebrate bait harvest. Transactions of the American Fisheries Society 122:46-62. Roseman, E.F., J.S. Schaeffer, E. Bright, and D.G. Fielder. 2014. Angler-caught piscivores diets reflect fish community changes in Lake Huron. Transactions of the American Fisheries Society in press. Roth, B.M., T.R. Hrabik, C.T. Solomon, N. Mercado-Silva, J.F. Kitchell. 2010. A simulation of food-web interactions leading to rainbow smelt Osmerus mordax dominance in Sparkling Lake, Wisconsin. Journal of Fish Biology 77:1379-1405. Royle, J.A., and R.M. Dorazio. 2008. Hierarchical modeling and inference in ecology: the analysis of data from populations, metapopulations, and communities. San Diego, Academic. Schick, R.S., S.R. Loarie, F. Colchero, B.D. Best, A. Boustany, D.A. Conde, P.N. Halpin, L.N. Joppa, C.M. McClellan, and J.S. Clark. 2008. Understanding movement data and movement processes: current and emerging directions. Ecology Letters 11:1338-1350. Schmalz, P.J., A.H. Fayram, D.A. Isermann, S.P. Newman, and C.J. Edwards. 2011. Harvest and exploitation. Pages 375-397 in B.A. Barton, editor. Biology, management, and culture of walleye and sauger. American Fisheries Society, Bethesda, Maryland. Schnute, J.T. 1994. A general framework for developing sequential fisheries models. Canadian Journal of Fisheries and Aquatic Sciences 51:1676-1688. 111 Schueller, A.M., M.J. Hansen, and S.P. Newman. 2008. Modeling the sustainability of walleye populations in northern Wisconsin lakes. North American Journal of Fisheries Management 28:1916-1927. Schueller, A.M., A.H. Fayram, and M.J. Hansen. 2012. Simulated equilibrium walleye production density under static and dynamic recreational angling effort. North American Journal of Fisheries Management 32:894-904. Schwarz, C.J., J.F. Schweigert, and A.N. Arnason. 1993. Estimating migration rates using tagrecovery data. Biometrics 49:177-193. Sethi, G., C. Costello, A. Fisher, M. Hanemann, and L. Karp. 2005. Fishery management under multiple uncertainty. Journal of Environmental Economics and Management 50:300-318. Spiegelhalter, D.J., N.G. Best, B.P. Carlin, and A. van der Linde. 2002. Bayesian measures of model complexity and fit. J. R. Stat. Soc. B, 64:583-639. Staggs, M. D., R. C. Moody, M. J. Hansen, and M. H. Hoff. 1990. Spearing and sport angling for walleye in Wisconsin's ceded territory. Wisconsin Department of Natural Resources, Bureau of Fisheries Management, Administrative Report 31, Madison. Steinhart, G.B., R.A. Stein, and E.A. Marschall. 2004. High growth rate of young-of-the-year smallmouth bass in Lake Erie: a result of the round goby invasion? Journal of Great Lakes Research 30:381-389. Strakosh, T.R., and C.C. Krueger. 2005. Behavior of post emergent lake trout fry in the presence of the alewife, a non-native predator. Journal of Great Lakes Research. 31:296-305. Thomas, L., S.T. Buckland, K.B. Newman, and J. Harwood. 2005. A unified framework for modeling wildlife population dynamics. Australian and New Zealand Journal of Statistics 47:19-34. VanDenAvyle, M. J. and R. S. Hayward. 1999. Dynamics of exploited fish populations. Pages 127-166 in C. Kohler and W. A. Hubert, editors. Inland fisheries management in North America 2nd edition. American Fisheries Society, Bethesda, Maryland. Vander Zanden, M.J. and J.D. Olden. 2008. A management framework for preventing the secondary spread of aquatic invasive species. Canadian Journal of Fisheries and Aquatic Sciences 65: 1512-22. Vander Zanden, M.J., and J.B. Rasmussen. 1999. Stable isotope evidence for the food web consequences of species invasions in lakes. Nature 401:464-467. Vander Zanden, M. J. and J. B. Rasmussen. 2001. Variation in δ15N and δ13C trophic fractionation: implications for aquatic food web studies. Limnology and Oceanography 46:2061-2066. 112 Vandergoot, C.S., and T.O. Brenden, T.O. 2014. Estimation of tag shedding and reporting rates for Lake Erie jaw-tagged walleye. Transactions of the American Fisheries Society 143:188-204. Vandergoot, C.S., T.O. Brenden, M.V. Thomas, D.W. Einhouse, H.A. Cook, and M.W. Turner. 2012. Estimation of tag shedding and reporting rates for Lake Erie jaw-tagged walleye. North American Journal of Fisheries Management 32:211-223. Vanderploeg, H.A., T.F. Nalepa, D.J. Jude, E.L. Mills, K.T. Holeck, J.R. Liebig, I.A. Grigorovich, and H. Ojaveer. 2002. Dispersal and emerging ecological impacts of PontoCaspian species in the Laurentian Great Lakes. Canadian Journal of Fisheries and Aquatic Science 59:1209-1228. Walters, C., and S.J.D. Martell. 2002. Stock assessment needs for sustainable fisheries management. Bulletin of Marine Science 70:629-638. Warburton, K., S. Retif, and D. Hume. 1998. Generalists as sequential specialists: diets and prey switching in juvenile Silver Perch. Environmental Biology of Fishes 51:445-454. Weeks, J. G. and M.J. Hansen. 2009. Walleye and muskellunge movement in the Manitowish Chain of Lakes, Vilas County, Wisconsin. North American Journal of Fisheries Management 29:791-804. Whitlock, R. and M. McAllister. 2009. A Bayesian mark-recapture model for multiple-recapture data in a catch-and-release fishery. Canadian Journal of Fisheries and Aquatic Science 66:1554-1568. Ying, Y., Y. Chen, L. Lin, and T. Gao. 2011. Risks of ignoring fish population spatial structure in fisheries management. Canadian Journal of Fisheries and Aquatic Sciences 68:21012120. 113