MULTISCALE INVESTIGATIONS INTO THE THERMAL HABITAT USE AND CONSERVATION OF ARCTIC GRAYLING (THYMALLUS ARCTICUS) IN THE PARSNIP RIVER WATERSHED, CANADA by Joseph Bottoms B.Sc., Minnesota State University, Mankato, 2013 DISSERTATION SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY IN NATURAL RESOURCES AND ENVIRONMENTAL STUDIES UNIVERSITY OF NORTHERN BRITISH COLUMBIA September 2024 © Joseph Bottoms, 2024 Abstract The drivers of the abundance and distribution of riverine species are often multiscale; they depend on both the local heterogeneity of habitats within the river and the context of the larger landscapes through which they flow. Methodological constraints can limit the ability of river ecologists to conduct studies over large scales, often producing results over fine scales that must then be extrapolated over unsampled units. Multiscale sampling approaches have been increasing in riverine studies, but the applications of their findings to conservation management – which operates over its own scales – are not always clear. In this dissertation, I investigate the use of physical and thermal habitats by Arctic grayling (Thymallus arcticus) at the reach, river, and watershed scales. Through experiments and observational studies, I investigated how this species used habitats in a degraded watershed. As a component of my research, I developed new statistical models that were parameterized with data from acoustic telemetry and applied new approaches using drone technology. Results from those models will allow fisheries ecologists to better understand the distribution of fish in river networks. I found strong associations between the distribution of Arctic grayling and pool habitats at the reach and river scales, and nonlinear relationships with temperature at the reach (11.0–16.0 °C) and watershed scales (11.1–17.1 °C). Unexpectedly, I found that the combination of upstream distance and pool habitats were a stronger predictor of Arctic grayling abundance than temperature at the river scale. This was explained by the fact that upstream distance accounted in part for the effects of temperature but also likely explained variation in other important predictors of Arctic grayling abundance (e.g. forage density, site fidelity, and territoriality). I quantified a thermal preference range (10.1–13.0 °C) for Arctic grayling in the laboratory and related this metric to the in-situ thermal habitat use of tagged but free-ranging individuals to determine the effectiveness of behavioural thermoregulation as a strategy for maintaining body temperature. I found that Arctic grayling used behavioural thermoregulation to effectively maintain their body temperatures during their summer feeding window, that heat transfer in Arctic grayling was slightly more efficient when cooling than warming, and that heat transfer was more rapid in smaller fish. I found novel evidence that Arctic grayling may be single-direction thermoregulators that will invest energy into cooling but not into warming. Contrary to my expectations, I found that the effectiveness of behavioural thermoregulation was more strongly related to thermal heterogeneity in time (along diel and seasonal axes) than in space. My findings related to the thermal ecology of Arctic grayling emphasized the continued conservation of coldwater-producing habitats (i.e. headwaters). I also found through my reach and river-scale studies that additional conservation of highly structured instream habitats that support thermal refugia and feeding opportunities across local scales would support better outcomes for the conservation of this population. I identified potential future risks, namely territoriality and dominance hierarchies that may keep Arctic grayling from leaving unfavorable thermal habitats during heatwaves. 3 Preface All animal capture and handling in this dissertation was done in accordance with BC Ministry of Forests scientific fish collection permits PG18-356580, PG19-523435, PH20-606121, PG21622265, and PG22-738069) and the University of Northern British Columbia’s Animal Care and Use Committee permits (ACUC protocols 2018-06 and 2021-05). Chapter 2 was developed as a manuscript to be submitted for publication with coauthors Marie Auger-Méthé, Bryce O’Connor, Michael Power, David A. Patterson, J. Mark Shrimpton, Steven J. Cooke, and Eduardo G. Martins. Chapters 3 and 4 were developed as an integrated technical report prepared for the Fish and Wildlife Compensation Program and were separated into standalone chapters for this dissertation. Both chapters are currently in development as their own manuscripts and follow the same format as chapter 2. These chapters will be submitted for publication with coauthors Eduardo G. Martins, Marie Auger- Méthé, Bryce O’Connor, John Hagen, Chris More O’Ferrall, and Avery Dextrase. 4 Table of Contents Abstract ........................................................................................................................................... 2 Preface............................................................................................................................................. 4 Table of Contents ............................................................................................................................ 5 List of figures .................................................................................................................................. 9 List of tables .................................................................................................................................. 15 Acknowledgements ....................................................................................................................... 19 Chapter 1 - Introduction ................................................................................................................ 23 Chapter 2 - Linking fish acoustic telemetry data to spatial covariates in river networks with spatial capture-recapture models .............................................................................................................. 34 Abstract ..................................................................................................................................... 34 Introduction ............................................................................................................................... 35 Methods ..................................................................................................................................... 38 Results ....................................................................................................................................... 46 Discussion ................................................................................................................................. 47 Chapter 2 Tables........................................................................................................................ 54 Chapter 2 Figures ...................................................................................................................... 57 Chapter 3 - The distribution of Arctic grayling (Thymallus arcticus) is related to temperature at the reach scale, but not at the river scale ...................................................................................... 64 Abstract ..................................................................................................................................... 64 5 Introduction ............................................................................................................................... 65 Methods ..................................................................................................................................... 69 Study Area ............................................................................................................................. 69 Snorkel surveys...................................................................................................................... 71 Drone surveys ........................................................................................................................ 73 Temperature covariate ........................................................................................................... 74 Other covariates ..................................................................................................................... 76 Modeling ................................................................................................................................ 77 Results ....................................................................................................................................... 79 Reach-scale results ................................................................................................................ 79 River-scale results.................................................................................................................. 81 Discussion ................................................................................................................................. 82 Reach-scale results ................................................................................................................ 83 River-scale results.................................................................................................................. 84 Chapter 3 Tables........................................................................................................................ 91 Chapter 3 Figures ...................................................................................................................... 94 Chapter 4 - Arctic grayling (Thymallus arcticus) show diel and seasonal patterns of behavioural thermoregulation ......................................................................................................................... 106 Abstract ................................................................................................................................... 106 Introduction ............................................................................................................................. 107 6 Methods ................................................................................................................................... 111 Study Area ............................................................................................................................ 111 Fish capture...........................................................................................................................112 Shuttlebox experiments ........................................................................................................... 113 Thermal preference experiments ..........................................................................................114 Heat transfer experiments .....................................................................................................114 Biologging with radiotelemetry .............................................................................................. 116 Quantification of thermal habitat ............................................................................................ 117 Data analyses ........................................................................................................................... 118 Heat transfer coefficient .......................................................................................................118 Predicting body temperatures of tagged fish ........................................................................119 Behavioural thermoregulation ............................................................................................. 121 Results ..................................................................................................................................... 124 Discussion ............................................................................................................................... 127 Chapter 4 Tables...................................................................................................................... 139 Chapter 4 Figures .................................................................................................................... 144 Chapter 5 - Discussion ................................................................................................................ 150 Methodological developments ................................................................................................ 151 Ecological insights and applications to conservation .............................................................. 155 References ................................................................................................................................... 165 7 Supplemental Information .......................................................................................................... 183 Appendix A. SCR Modeling tutorial .......................................................................................... 216 Introduction ............................................................................................................................. 217 Read and inspect prepared data files ....................................................................................... 218 Read traps files. ................................................................................................................... 218 Read detection data .............................................................................................................. 219 Read covariate data .............................................................................................................. 220 Standardize covariates for analysis ......................................................................................... 221 Create masks ........................................................................................................................... 221 Sampling effort ........................................................................................................................ 222 Define capture history ............................................................................................................. 224 Fit candidate SECR models to data ......................................................................................... 225 AIC selection of the best model .............................................................................................. 226 Predict the temperature covariate from the top model ............................................................ 227 Plotting model spatial outputs ................................................................................................. 229 Data prep .............................................................................................................................. 229 Create plots .......................................................................................................................... 230 Save outputs for visualization in GIS .................................................................................. 232 Appendix B. Expanded radio tag inspection plots ...................................................................... 234 8 List of figures Figure 1.1. A male (left) and female (right) Arctic grayling (Thymallus arcticus). Photo by Pat Clayton courtesy the Center for Biological Diversity, available through public domain. ............ 28 Figure 1.2. The Williston Reservoir watershed, with the Parsnip River watershed highlighted in orange. Work in this study occurred in the Parsnip River mainstem and its four major tributaries: the Anzac, Table, Hominka, and Missinka Rivers........................................................................ 30 Figure 2.1. The Parsnip River watershed in north-central British Columbia, Canada, was discretized into a state-space mask of 366 pixels at a resolution of 1 Rkm/pixel. The state-space mask includes the Parsnip River mainstem and four major tributaries included in this study: the Anzac, Table, Hominka, and Missinka Rivers. The watershed flows northwest through the Rocky Mountain Trench into the Williston Reservoir. Acoustic receiver deployments are indicated by year; overlapping points indicate replacements of compromised or lost acoustic receivers. ....... 58 Figure 2.2. Model-averaged predictions of the relative density of tagged Arctic grayling activity centres with relation to the temperature covariate across all three summer seasons (panel A) and the bull trout covariate BT in the spring and summer (panel B). Relative density of tagged Arctic grayling peaked at 11.3 °C (25th and 75th quantiles 11.1 – 17.1 °C). The shaded region indicates the ± 1 standard error of model predictions. Tagged Arctic grayling relative density response to tagged bull trout relative density was positive in both the winter and summer, with a stronger response in the winter. Areas where the distributions of standard error width is narrow have few data points of activity centre density estimates at those temperatures in 2019 and 2020. ............ 59 9 Figure 2.3. Annual summer distributions of tagged Arctic grayling activity centres estimated by the SCR model across the Parsnip River watershed. Activity centres are given by white markers over the relative tagged population density across the state-space in blue. Coloured circles indicate the receiver array and which acoustic receivers were active over the detection window (Jul. 1 – Sep. 15 each year). Summer activity centre distributions were wide ranging within the Parsnip River and its four main tributaries with mainstem Parsnip River activity centres being closely associated with the Table River confluence.................................................................................. 60 Figure 2.4. Annual spring distributions of tagged Arctic grayling activity centres estimated by the SCR model across the Parsnip watershed. Activity centres are given by white markers over the relative tagged population density across the state-space in blue. Coloured circles indicate the receiver array and which acoustic receivers were active over the detection window (May 1 – Jun. 1 each year). Spring activity centre distributions occurred largely in the Parsnip River mainstem and the lower to middle reaches of the tributaries. ....................................................................... 61 Figure 2.5. Annual winter distributions of tagged Arctic grayling activity centres estimated by the SCR model across the Parsnip watershed. Activity centres are given by white markers over the relative tagged population density across the state-space in blue. Coloured circles indicate the receiver array and which acoustic receivers were active over the detection window (Nov. 1 – Mar. 31 each year). Winter activity centre distributions were concentrated to the Parsnip River mainstem and the lower reaches of the Anzac River. One activity centre appearing in the upper Missinka River in 2021 is suspected to be a case of mortality as the tag ID does not reappear in the subsequent spring or summer models for 2021. ........................................................................... 62 Figure 2.6. Model-averaged spatially-explicit predictions of the density of Arctic grayling activity centres in the Parsnip River watershed based on the mean temperature at each location in the state 10 space in the 2019 (a warm year) and 2021 (a cool year) summer seasons. Distribution is predicted to be higher in the tributaries and upper Parsnip River where summer temperatures are cooler than the Parsnip River mainstem. Predicted densities show wide Arctic grayling distributions during the cool year and become patchier and more discrete during the warm year. .............................. 63 Figure 3.1. Orthomosaic image of the Anzac River from river kilometers (Rkm) 0 - 51. This figure is broken into three segments for the lower (left), middle (center), and upper Anzac River (right). River segments have been rotated per the corresponding compass orientations to fit in one image. Purple regions highlight the two nested study reaches in this study............................................. 95 Figure 3.2. Orthomosaic image of reach “A”, the lower of the two reaches from the reach-scale study from Rkm 37 - 39. It features two small tributary inputs, the lower of which enters the Anzac at temperatures warmer than mainstem temperatures after flowing through regions of clearcut land use (Callout A). Habitats in reach A were varied, including cobble-boulder pools along a bedrock shelf, a long shallow run over bedrock substrates, a prominent deep corner pool (Callout B), river braids around instream islands, and multiple series of riffle-pool interfaces (Callout C). A long section of shallow riffles and micro-pools spans through the lower parts of the reach (Callout D) before opening to small shallow pools and tailing out of the study area. Subreaches from the snorkel surveys (e.g. AR01) are labeled by the black outlines with white text. Imagery has been slightly rotated per the compass rose from true north to fit on the page. ..................................... 96 Figure 3.3. Orthomosaic image of reach “B”, the upper of the two reaches, was a 1.1 km reach defined from Rkm 46.5 - 47.6. It featured a major tributary input which supported some Arctic grayling use in the lower 100 m (Callout C), a major bedrock chute and series of rapids (Callouts 11 A, B), series of cobble riffles and pools (Callout D), a prominent plunge pool below the chute (Panel A), and a prominent feeding pool below the tributary confluence (Callout C). Subreaches from the snorkel surveys (e.g. BR01) are labeled by the black outlines with white text. ............ 98 Figure 3.4. The five size classes of pool habitats defined at the river scale. Size 0: Very small, but larger than micropools which cannot be clearly identified at the imagery resolution. Possibly ephemeral features caused by temporary river features such as isolated logs or exposed large boulders creating collecting pools in their wake. Took up < 10% of the overall reach width; Size 1: Small pools, approximately equal in length and width, taking up 10 - 25% of the overall reach width; Size 2: Pools which took up 25 - 100% of the overall reach width but were less than two reach widths in length and/or contained inclusions of other habitat types; Size 3: Large pools which took up the entire width of the reach and had a length of two reach widths or more; Size 4: Notably large pools or pool complexes which covered the entire reach for three or more reach widths or which showed major widening of the reach when compared to adjacent channels. ................... 100 Figure 3.5. Model averaged predictions of the effects of temperature (Panel A) and pool size (Panel B) on Arctic grayling abundance at the reach scale. An interaction effect between these two predictors was found, with a peak response in abundance at temperatures of 12.5 °C in large pool habitats, 13.5 °C in medium pool habitats, and 15.3 °C in small pool habitats (Panel B). The response in abundance to pool areas was highest at cool temperatures, weakly increased with pool area at medium temperatures, and weakly declined with pool area at high temperatures. The three panels of each series of plots represent the 0.25, 0.50, and 0.75 quantiles of temperatures and pool sizes, respectively. Plot points indicate the observed data in each modeled combination. Shaded regions represent the 95% CI for the predictions........................................................................ 102 12 Figure 3.6. Model averaged predictions of Arctic grayling abundance as a function of river kilometer from the reach-scale study. Reach B showed a stronger response with higher Arctic grayling counts than Reach A. Plot points indicate the observed data from the snorkel surveys. ..................................................................................................................................................... 103 Figure 3.7. Model-averaged predictions of Arctic grayling abundance as a function of river kilometer (Rkm; Panel A), mean pool size per Rkm (Panel B), and pool density per Rkm (Panel C) from the river-scale study. The strongest predictor was Rkm with mean pool size showing a similarly strong effect, whereas the effect of pool density was small. ....................................... 104 Figure 3.8. Model averaged predictions of Arctic grayling abundance in the Anzac River at a 1 Rkm resolution. The two models used for predictions included terms for upstream distance (Rkm), mean pool size per Rkm (top and second ranked models), and pool density per Rkm (second ranked model only). Trends in predicted Arctic grayling abundance were low in the lower Anzac River with a few patches of higher densities around Rkms 8 and 17 and increased with upstream distance. Abundance predictions were highest in the upper river with the highest predictions occurring around Rkm 45 and decreasing slightly into the highest reaches. .............................................. 105 Figure 4.1. Distribution of temperature loggers and radio receivers used in this study. ............ 144 Figure 4.2. Predictions of heat transfer in Arctic grayling along a body weight gradient as exposed to a heating (12 – 15 °C) and cooling (15 – 12 °C) temperature shift. Smaller fish were quicker to heat or cool than larger fish, but larger fish equilibrated at a slightly higher temperatures due to 13 higher metabolic heat production. The rate at which body temperatures equilibrated was similar, with modeling indicating slightly more efficient cooling than heating (Table 3). ..................... 145 Figure 4.3. Regression analysis of the index of accuracy of thermoregulation, ܾ݀, against the index of thermal quality of the habitat, ݀݁. The dashed line represents the line of thermoconformity (slope = 1). Points above this line indicate ineffective thermoregulation and points below this line indicate effective thermoregulation. The regression line of the points falls below the line of thermoconformity and indicates thermoregulatory behaviours. ................................................. 146 Figure 4.4. Regression analysis of the degree of hourly thermal exploitation ‫ݔܧ‬, or the proportion of time Arctic grayling occupied their thermal preference range (ܶܵ‫ )ܶܧ‬when it was available in their environments. ..................................................................................................................... 147 Figure 4.5. An example plot from Arctic grayling no. 15 showing (A) body temperature in blue against median water temperatures in red with the black horizontal lines denoting ܴܶܲ‫( ܨܧ‬solid) and the upper and lower bounds of ܶܵ‫( ܶܧ‬dashed) against study hour, (B) the accuracy of body temperatures with respect to ܶܵ‫ܾ݀ ܶܧ‬, (C) index ‫ ߂ܧ‬as positive or negative departures from zero (thermoconformity) with positive values representing thermoregulatory behaviours and negative values denoting avoidance of preferred temperatures and day or night period by color (a general trend from day to night thermoregulation can be seen starting around hour 850 on September 4), and (D) the index of thermal exploitation (‫ )ܺܧ‬against study hour by day and night. Plots for all fish with sufficient data can be found in Appendix A. ............................................................... 148 Figure 4.6. Model-averaged prediction plots for the GAMM analysis of the thermoregulatory effectiveness (index ‫ܧ‬Δ) against (A) mean hourly activity, (B) body condition, (C) the patchiness 14 of the thermal habitat, and (D) hour of the day by season. The dotted horizontal line represents the line of thermoconformity, and the shaded areas bounded by the dashed lines represent the standard error of the predictions. ............................................................................................................... 149 List of tables Table 2.1. Annual counts of Arctic grayling Thymallus arcticus (GR) and bull trout Salvelinus confluentus (BT) acoustically tagged in the Parsnip River and its four major tributaries. The species codes displayed correspond to the codes used by the Ministry of Environment for British Columbia. ...................................................................................................................................... 54 Table 2.2. Annual deployments and active coverage extent of the acoustic receiver array by river. High receiver loss occurred during a high discharge freshet in 2020. See Figure 2.1 for a visual description of these data. Values in parenthesis represent the receiver deployments (+) and losses (-) during each year. ...................................................................................................................... 55 Table 2.3. AIC statistics for candidate Arctic grayling multisession SCR models for each season. Presented are the log likelihood, the standard AIC score and the AICc score for small sample sizes, the difference in AICc between each model from the top model (ΔAICc), and the respective model weights (AICcwt). Model terms are the relative density response (D), temperature and squared temperature, and the bull trout covariate BT. ............................................................................... 56 15 Table 3.1. Table of each drone and snorkel survey conducted in the reach-scale study within the Anzac River, as well as the calibration values between the FLIR sensor and HOBO temperature loggers at the time of the flight. .................................................................................................... 91 Table 3.2. AICc table for the reach-scale modeling study, where GR is the count of Arctic grayling in each 100-m subreach, ‫ ݐ‬+ ‫ݐ‬2 is the polynomial temperature term, ℎ. ‫ ݔ‬is the area in ݉2 of each habitat type ℎ. ‫( ݔ‬pool, riffle or run) within each subreach, ܾ is the count of bull trout during each survey, ‫ ݓ‬is a binary presence/absence parameter of mountain whitefish occurrence in each subreach during each survey. All models include a random effect of subreach on the intercept. Presented are the number of parameters in the model (K), the log likelihood (logLik), the AIC score with small sample size correction (AICc), the difference in AICc from the top model (dAICc), and the model weights (AICcwt). Models which are bolded in ‫ ܶܧܵܯ‬2 were included in the 95% confidence set. ................................................................................................................ 92 Table 3.3. AICc table for the river-scale modeling study, where ݈ is the offset logged length of each snorkel survey in km to account for the variable length (1.2 – 5.4 Rkm) of the snorkel surveys (݈ is included in all models), ‫ ݑ‬is the upstream distance from the confluence of the Anzac River with the Parsnip River in km, ‫ ݐ‬+ ‫ݐ‬2 is the polynomial temperature term, ‫ ݏ‬is the mean size of all pools within each Rkm, and ݀ is the density of pools within each Rkm. Presented are the number of parameters in the model (K), the log likelihood (logLik), the AIC score with small sample size correction (AICc), the difference in AICc from the top model (dAICc), and the model weights (AICcwt). Models in bold were included in the 95% confidence set. .......................................... 93 16 Table 4.1. Experimental design for the shuttlebox heat transfer experiments. Each of the five experimental trials (shifts) exposed an individual Arctic grayling to a warm (+) or cool (-) temperature differential. The initial direction of the temperature shift pattern was determined by a randomized start in either the cool tank or the warm tank. As trials were conducted, the magnitude of the temperature differential gradually decreased from a starting differential of 2.0 °C to a final differential of 1.0 °C. .................................................................................................................. 139 Table 4.2. Summary table for the 10 Arctic grayling used in the shuttlebox thermal preference experiments. The thermal preference (ܴܶܲ‫ )ܨܧ‬is given as the median (q50) of occupied temperatures (ܱܶ‫ )ܥܥ‬during the twenty hours of data collection and the thermal preference range (ܶܵ‫ )ܶܧ‬is given as the temperatures in the 25th (q25) and 75th (q75) quantiles. ........................ 140 Table 4.3. AICc statistics for the non-linear mixed effects models of the heat transfer coefficient ݇ against centred body weight (Wc) and whether the body temperature of the fish was warming or cooling relative to ambient water temperatures (trend [warm]). Presented are the log likelihood, the number of parameters in the model (K), the AICc score, the difference in AICc between each model from the top model (ΔAICc), and the respective model weights (AICcwt), and the 95% confidence set (highlighted in bold). A random effect was included in the model for each fish used in the experiment......................................................................................................................... 141 Table 4.4. AICc statistics for the candidate GAMMs modeling the index of thermoregulatory effectiveness ‫ ߂ܧ‬against hour of the day (hour), body condition (bodycon), activity (act), season (before and after September 4), and patchiness (patch) of the thermal habitat. Presented are the log likelihood (logLik), the number of parameters in the model (K), the AICc score for small sample sizes, the difference in AICc between each model from the top model (ΔAICc), the respective 17 model weights (AICcwt), and the cumulative model weights (cAICcwt) for selecting the 95% confidence set for the best model (highlighted in bold) ............................................................. 142 Table 5.1. Insights into Arctic grayling ecology found in this dissertation and their implications into the conservation of the Parsnip River watershed population. The umbrella term ‘Protection’ is used to describe conservation implications, but this can encompass habitat restoration, enhancement, land procurement or protection actions. .............................................................. 155 18 Acknowledgements An enormous number of people helped me to finish this dissertation over the past five years. To Dr. Eduardo Martins, thanks for taking me on as a student, encouraging me to take the chance on rolling up to the PhD, and providing expert guidance through this process from start to finish. Your mentorship and friendship (through all my varied levels of grumpiness) has been deeply appreciated, from catching the little things in my code that I was too close to see, to providing the research tools and direction for me to make this project my own, to the occasional neighborly lending of power tools and backup drysuits. Your constant support from start, finish, and beyond grad school have been life-changing and I will be forever grateful. To Dr. Marie Auger-Méthé, thank you for jumping in with both feet to co-supervise me through this process and throughout it changing my understanding of what it means to do next-level science. Your patient reviews and workshopping of my early drafts where I attempted to use Cadillac statistics when I was still riding with training wheels have me eternally in your debt. Your kindness and attention when your own plate was overflowing have been so meaningful, and the gently delivered and honest reality checks have kept me from making huge mistakes when loading my own plate. Thank you to my committee members Drs. Mark Shrimpton and Chris Johnson. Dr. Shrimpton’s knowledge of this system both past and present has been instrumental in shaping the direction of this work, and Dr. Johnson’s advice that research does not work in widgets has been thoroughly internalized and will not be forgotten. Thank you to Dr. Pam Wright who held a brief tenure on my 19 committee and provided a crash course in effective conservation planning before moving on to a well-deserved retirement. To my lab cohort and friends, the coffee breaks and group chats have been a fantastic source of sanity and lighthearted commiseration over the years. Individually, you’ve all helped me get through this in ways I feel like I have never been able to fully repay. Thanks to Annika Putt for the continual advice, asking the hard questions during practice sessions, and for bringing forwardthinking into our lab extracurriculars. Thanks to Avery Dextrase for being the glue holding everyone together, always supporting your labmates in the most unfailing way, for teaching me how to radio tag, troubleshoot exploding chillers, avoid cultural faux-pas, and for offering a driveway to keep my field camper safe for the winter. Thanks to Bryce O’Connor for the advice on fly fishing when it was still new to me, without which the latter half of the project would have failed completely, and for the advice about working in the Anzac, which saved us untold hours of walking. Thanks to Liz Hirsch, who is an absolute riot to work with, who helped me to make sense of the Frankenstein’s monster known as the shuttlebox, and whose dark poetry about driving the harrowing Crocker FSR still hangs on my fridge. Thanks to Meaghan Rupprecht, who has been my geospatial kindred spirit and was always a wonderful sounding board for troubleshooting the fantastic yet frustrating process of making maps in R, who helped me pick up the pieces of a mild transport disaster along the Hart highway, and who sacrificed some personal dry flies to the success of the project in the eleventh hour (though maybe they were Nigel’s? Thank you Nigel). Thanks to Devon Smith, a steady friend who started as a field tech on this project and is now close to finishing his own thesis, whose sense of organization and direction kept me alive despite my navigational impairments as we travelled Europe during a telemetry conference, and whose dedication in the field was seemingly never phased by overlanding deep-cycle batteries. Thanks to 20 Dan Larson, who as a fellow Minnesotan was an immediately friendly face and a compelling force of laughter when I first started, and who unfortunately left us too soon. Thanks to my newer labmates, too, who came along after I finished the work but whose friendship has been valuable and I look forward to seeing what you do in the coming years; Lucas Moura, Eliseu Gomes, Abby Oviatt, Melody Mah, and John Gray. Thanks to the members of my remote lab-away-from-home under Marie’s supervision in the Statistical Ecology Research Group for always being friendly and welcoming, even if we’ve never had the chance to meet face-to-face. To the field technicians who helped get this work done in the field, I promise it was never my intention to overwork you. Thanks to Chris More O’Ferrall, who was an instant lifelong friend whose experience in the bird world got my head out of the water and kept the drone project alive, whose field experience let me pivot the project a hundred times a day without skipping a beat, whose easy humor let the long, wet days carry on with laughter, and who could seemingly patch a tire with anything. Thanks to Ian Clevenger, who could always make the long days jetboating the Hominka River searching for acoustic receivers fun, and who understood that the Bear Lake gas station is truly the hanging gardens of our time. Thanks to Danny Scurfield, who could operate a quad like only a proper Alberta boy could, and whose thoughtful perspectives on natural resource issues was always refreshing. Thanks to Grayson Vanderbyl, whose early prototyping of remote shuttlebox work and utter dedication to camp living set the groundwork for future work during my first field season. Thank you to those in my life who supported me from outside of the academic setting. Thanks to Kari Van Ruskenveld, who offered steadfast support through all the dissertation highs, lows, and concussions, who always seemed to have the perfectly complementary skillset for the things I was missing, who made sure I was fed and watered during all the late-night pushes, who with no 21 hesitation got away and explored B.C. with me when I needed a break, who shuttled field techs to and from the Anzac and filled in when I had nobody to go with me, and who always reminded me that a day of failed progress was still progress in its own way. Thank you to my parents Bart and Lisa, who stayed excited about the work I was doing and helped develop, build, and deliver my truck camper which kept me warm, comfortable, and sane during the extensive fieldwork. Thanks to Evan and Mitch Ogaard and Danielle Grunzke for carrying on with me even when I dropped off the radar for months at a time. Thanks to Andrea Aspelund, who was beyond excited about this opportunity for me and who we unfortunately lost before it saw its conclusion. Thank you to those who provided the funding for this project; the Fish and Wildlife Compensation Program, the Canadian Foundation for Innovation, Natural Sciences and Engineering Research Council of Canada, and the UNBC Research Project Award. This support was the only way that any of this was possible. Thank you to the First Nations whose traditional territories this work was conducted on: Treaty 8 signatory Nations Prophet River, Saulteau, and West Moberly First Nations, the McLeod Lake Tse’Khene, and Lheidli T’enneh First Nation on whose traditional and unceded territory I live and study. To the countless others who helped by offering technical advice, help in the field, or just overall friendly support, thank you. Any attempt to compile an exhaustive list of everyone who helped me along the way is almost certain to be incomplete, so I apologize if anyone is missing. With that said, thank you to Matthew Blackburn, Dr. Murray Efford, Dr. Chris Sutherland, John Hagen, Mike Stamford, Dr. Ellen Petticrew, Tom Balfour, Tom Willms, Paige Wilson-Breault, David Breault, Brittney Reichert, Amber Pallister, Sydney Angielski, Mackenzie Howse, Carly Walters, Emily Mason, Julian Brown, Rachelle Foubert, Julien Wauthier, Jerome Egger (gone too soon), Dr. Nikolaus Ganter, Ian Spendlow, Riognach Steiner, Julian Napoleon, Cale Babey, Dr. Daniel Sims, 22 Chris Morgan, Luc Turcotte, Andy Rosenberger, Jeannine Randall, Matt McLean, Dr. Joe Shea, Amanda Carriere, Jaqueline Middleton, Rebecca Lahti, Ainsley Davison, Jesi Lauzon, Joel KhooEn, Alex Bevington, Seb Dalgarno, Dr. Dezene Huber, and Dr. Danie Erasmus. Oh, and of course, my therapists. Chapter 1 - Introduction Understanding how fish access and exploit suitable habitats within degraded landscapes is important to conservation planning in the Anthropocene (Veech 2021). As climate change forces more frequent and more severe climatic extremes, effective conservation management will require a further understanding of how animals modify their behaviours to offset the challenges set forth by their environments. While research objectives are defined based on our gaps in knowledge, addressing these objectives across spatiotemporal scales can often be constrained by the tools and techniques available to scientists for studying animal populations. These can be related to a mismatch in data types and existing statistical approaches, scalability issues in data collection or processing timelines, or changing abiotic conditions within the study systems. Problems presented by scale have been long examined in ecology (Wiens 1989; Levin 1992). Studies at small spatial scales can identify the drivers of animal habitat selection but lack the context of how these drivers are connected across the greater landscape (Fausch et al. 2002). In contrast, studies at large spatial scales can be useful for coarse management purposes, but can miss discrete, but critical, habitat features (e.g. Fausch et al. 2002; Torgersen et al. 1999). For example, thermal refugia can drive animals towards discrete habitat patches that would not be accounted for when making predictions across coarse scales. Further, empirical studies in ecology generally occur over relatively small temporal scales, while emergent issues related to climate change 23 progress over decades (Schnieder 2001). At times, these constraints can limit the findings of empirical studies from being relevant to the scales at which conservation managers can affect change (Crossin et al. 2017). How animals adapt to a changing world is a complex interplay between the short-term behavioural plasticity of individuals (over minutes to weeks), the longterm evolutionary responses of populations (over generations; Kelly 2019; Penney et al. 2022), the heterogeneity of the environments in which they live (across small and large spatial scales; Sears et al. 2019; Sears and Angiletta 2016), and how this heterogeneity overlaps with the environmental preferences of individuals. While marine fishes have been observed moving to higher latitudes at an increasing frequency in response to climate change (Dahms and Killen 2023), the geographical constraints imposed on freshwater fishes by their habitats may limit this option across appropriate latitudinal scales (Busch et al. 2012). The importance of heterogeneous habitats is thus underscored in freshwater environments by the frequent inability of species to move northwards to escape climate impacts (excepting systems which flow north/south; Busch et al. 2012). In lakes, vertical migration within the water column can offset some climatic effects (Busch et al. 2012), though the turbulent mixing and often shallower depths of river habitats can reduce the availability of vertical clines in which to exploit this strategy (Torgersen et al. 1999; Dugdale 2016). Species occupying rivers may have the option to move to higher altitudes in systems where topography permits (Comte and Grenouillet 2013), though headwater reaches can be inaccessible due to steep gradients, narrow channel widths, and physical barriers (Isaak and Young 2023). Consequently, the heterogeneity of horizontal habitats may be one of the best options available for freshwater fishes in rivers to mitigate the effects of climate change. 24 Temperature has been characterized as the ‘master variable’ of freshwater fish ecology (Brett 1971). As with other ectotherms, temperature can act as a cue for migration (Elsner & Shrimpton 2019), spawning (Hubert et al. 1985; Shuter et al. 2012), and diel foraging patterns (Armstrong et al. 2013), dictate the times conducive to (and intensity of) activity (Gunderson & Leal 2016; Abram et al. 2017), and influence survival and reproductive success (Martins et al. 2011; Dahlke et al. 2020). Temperature distributions in rivers are driven by a suite of abiotic factors including atmospheric conditions (temperature, insolation, precipitation), discharge (volume, turbulence), streambed (hyporheic exchange, sediment properties), and topography (riparian land cover, aspect, slopes, surrounding landscapes; Caissie 2006; Fausch et al. 2002). The horizontal heterogeneity of temperature in rivers can vary across both small spatial scales (10s – 100s of m; Kurylyk et al. 2015; Dzara et al. 2019) and larger scales (headwaters to mouth). Temperature gradients can be driven by the longitudinal connectivity of fluvial habitats known as the River Continuum Concept (Vannote et al. 1980; Fausch et al. 2002). The response of freshwater fishes to climate-driven warming of riverine habitats will be related to both their physiological responses to warming and the heterogeneity of thermal habitats that enable thermoregulation (Huey and Slatkin 1976; Blouin-Demers and Nadeau 2015; Wolkovich et al. 2014). The physiology of freshwater fishes is mediated by temperature through their thermal optimum, or the range of temperatures (between the upper and lower pejus temperatures; Pörtner and Farrell 2008) at which the performance of metabolically driven activities (e.g. swimming speed, digestion rate, aerobic scope) are maximized (e.g. Eliason et al. 2011; Farrell et al. 2008; Anttila et al. 2013). Current thermal optima in salmonid populations are connected to historic (long-term) trends in their thermal environments (Farrell et al. 2008), and responses to temperatures that exceed these optima is both cumulative (on a daily scale) and proportionate to the magnitude and duration (over 25 a relatively fine scale) of the thermal stress (Farrell et al. 2008; Rezende et al. 2014). Freshwater fishes can recover from the cumulative effects of moderate thermal stress by spending adequate time within recovery temperatures. However, exposure to extreme temperatures that surpass physiologically rigid upper temperature thresholds, can trigger a cascade of physiological effects that can quickly become fatal (e.g. Pörtner 2002; Farrell et al. 2008; Martins et al. 2011). In real-world settings, the difference between an animal’s optimal thermal range, a metric defined relative to performance or its fundamental thermal niche, and its realized thermal niche, the actual temperatures occupied by free-ranging individuals (Huff et al. 2005; Ángeles-González et al. 2020) are driven by the thermal preference of an individual, or the range of temperatures selectively occupied by an individual within a free-choice environment (Ángeles-González et al. 2020). Freshwater fishes residing in heterogeneous thermal habitats can use behavioural thermoregulation (i.e. moving between warmer or cooler habitats) to maintain their internal body temperatures within their preferred thermal range (Haesemeyer 2020; Amat-Trigo 2022). Thermal preference ranges of individuals can differ from thermal optimums, but with some consequences to survival and reproduction. For example, theoretical models by Martin and Huey (2008) found that individuals occupying temperatures on the high end of thermal optimum curves had negative fitness outcomes relative to those occupying temperatures below curve peaks. Therefore, the thermal preference of freshwater fishes may be cooler than the optimum suggested by thermal performance curves (e.g. Larsson 2004). As such, thermal preference may be more relevant for conservation management than fundamental thermal niche with respect to characterizing available thermal habitats across the landscape (i.e. preference is optimal). The efficiency of thermoregulatory behaviours (i.e. how effectively an individual can maintain its body temperature) by individuals is driven by the relationship between their thermal preference 26 range, realized thermal niche, and the heterogeneity of thermal habitats (Sears et al. 2016). In practice, an animal may not always occupy its preferred thermal habitats, even when they are available to them. It is well understood in habitat ecology that an animal will occupy non-preferred habitats if the benefits of doing so outweigh the costs (though the costs can accrue to influence excursion length into non-preferred habitats; Veech 2021). In riverine habitats, favorable thermal habitats may be disjointed through space and time from areas with suitable foraging opportunities, flow rates and refugia, or predation and competition (Brennan et al. 2019). Using behavioural thermoregulation, freshwater fishes can balance the maintenance of body temperatures relative to these tradeoffs across fine temporal scales. While behavioural strategies may be suitable for offsetting the effects of climate change over the short term, behavioural thermoregulation may inadvertently suppress the adaptive capacity of populations over the long term (Huey et al. 2012; Kelly 2019; Penney et al. 2020). Behavioural thermoregulation may be further influenced by the presence of antagonistic species. For example, species may opt to use less suitable thermal habitats if it reduces their risk of predation (Webb and Whiting 2005). All these factors must be considered in the development of conservation plans (e.g. Blackman 2001; Stamford et al. 2017; Fisheries and Oceans Canada 2020) for freshwater fishes. For my dissertation, I demonstrate how to apply multiscale ecological insights to the conservation of fluvial Arctic grayling (Thymallus arcticus) in the Parsnip River watershed of north-central British Columbia, Canada. Arctic grayling are cold-water salmonids characterized by a prominent dorsal fin with vivid orange and blue coloration which is most pronounced in spawning males (McPhail 2007; Figure 1.1). They can follow three life history patterns: fluvial (most common), lacustrine, and adfluvial (migrations between streams and lakes; McPhail et al. 2007). Fluvial Arctic grayling undergo complex annual migrations between summer feeding, spring spawning, and overwintering 27 locations (McPhail 2007; Blackman et al. 2002). They are slow-growing and slow maturing, typically spawning for the first time between 4 – 7 years old and reaching sizes of 500 – 600 mm (McPhail 2007). They are unique among salmonids in that males arrive to defend spawning sites ahead of females (McPhail 2007). Figure 1.1. A male (left) and female (right) Arctic grayling (Thymallus arcticus). Photo by Pat Clayton, courtesy Center for Biological Diversity, available through public domain. Following the 1967 impoundment of the Williston Reservoir, the fluvial Parsnip River watershed population of Arctic grayling experienced significant habitat loss as flooding converted extensive areas of riverine habitats into deeper, slower, more lake-like habitats (Lashmar and Ptolemy 2002; Clarke et al. 2007; Stamford et al. 2017, Hagen and Gantner 2019). The effects of habitat loss were compounded by overfishing as expanding networks of resource roads increased anglers’ access to Arctic grayling habitats which were previously considered remote (Lashmar and Ptolemy 2002). By 1995, a harvest moratorium was imposed on the Williston watershed Arctic grayling following 28 their designation as a red-listed population (Northcote 1993; Lashmar and Ptolemy 2002). The use of Conservation Units in British Columbia was in its infancy in 1995, and as the program matured and criteria for genetic designation became established, the Williston watershed (Figure 1.2) Arctic grayling were reclassified as part of the greater Southern Beringean genetic lineage and moved to the yellow-list ca. 2002 (M. Stamford, Stamford Environmental; S. Pollard, Freshwater Fisheries Society of B.C. Personal Communications). Understanding how Arctic grayling use their thermal habitats in space and time is a critical data gap for the implementation of long-term conservation actions (Stamford et al. 2017). As they are now isolated from the rest of the Arctic watershed drainage at the southern periphery of their distribution in B.C. and separated from the rest of the Southern Beringean lineage by the Williston Reservoir, the effects of climate change on this population may be particularly acute over both behavioural and genetic time scales (Vatland et al. 2015; Troia et al. 2019). As the climate warms, Arctic grayling are likely to face summer energy deficits and in turn will have to allocate more resources towards behavioural thermoregulation or spend more time in thermal refugia. These short-term behaviours will have to be balanced against their energy requirements for growth and reproduction to ensure long-term success (Armstrong 2021). At both fine and coarse spatial scales, the distribution of habitats that support preferred temperatures will be important. 29 Figure 1.2. The Williston Reservoir watershed, with the Parsnip River watershed highlighted in orange. Work in this study occurred in the Parsnip River mainstem and its four major tributaries: the Anzac, Table, Hominka, and Missinka Rivers. 30 In this dissertation, I examine the habitat use of Arctic grayling across variable spatial and temporal scales. My objectives are to learn what thermal and physical habitat characteristics Arctic grayling rely on during their critical summer feeding window, how these are distributed within their seasonal habitats, and how they use behavioural strategies to interact with their thermal environments. In the second chapter, I refined established techniques into a novel way to link acoustic telemetry data to spatial covariates. I modified spatial capture-recapture models (Royle et al 2014; Efford 2023), which have been commonly applied in the terrestrial environment, to complex branching riverine systems. With the help of Dr. Murray Efford (University of Otago), I further modified these models to make them suitable for acoustic telemetry data. Then, I used the models to explore how the seasonal distribution of Arctic grayling is driven by both spatial variation in temperature and predation risk from bull trout (Salvelinus confluentus) at the watershed scale (10s – 100s of km). My third chapter characterized both thermal and physical habitat use of Arctic grayling in the Anzac River by combining counts of fish from snorkel surveys with available habitat distributions from drone mapping at both the river (1 – 52 km) and reach (100 m – 2 km) scales. My fourth chapter examined the thermoregulatory behaviour of Arctic grayling themselves. In it, I defined a thermal preference range for adult Arctic grayling during their summer feeding period in the Anzac River and the rate at which Arctic grayling body temperatures equilibrate after being exposed to a shift in ambient water temperatures. I used temperature-sensing radio telemetry to relate these metrics to the in-situ habitat use of free ranging individuals to calculate metrics of thermoregulatory effectiveness (Blouin-Demers and Weatherhead 2001; Blouin-Demers and Nadeau 2005). Together, these studies occured over three spatial scales (watershed, river, and reach) across three temporal extents (three years, one snapshot of physical habitats compared against three years of snorkel data, and four bi-weekly surveys, 31 respectively). I looked at habitat use at the population level in Chapters 2 and 3 (specifically, sampled population dynamics, which is of particular relevance to Chapter 2), and at the individual level in a study of thermal preference and habitat use in Chapter 4. Through these three chapters, I explored four general hypotheses (with more specific objectives defined in each chapter). As Arctic grayling are coldwater specialists, I hypothesized that summer water temperature will be a suitable predictor of their distribution across all three spatial scales examined in this dissertation. I hypothesized that to maintain their body temperatures around their thermal preference range, Arctic grayling will use behavioural thermoregulation, and that these strategies will be energetically costly to maintain during the hot summer months. As fluvial Arctic grayling are often associated with pool habitats (e.g. Blackman et al. 2002; McPhail 2007), I hypothesized that pool habitats will be important to explaining the variation in Arctic grayling distributions at the reach and river scales (though sampling this dynamic at the watershed scale would be prohibitive). I expected that at multiple spatial scales, the presence of potentially antagonistic species (bull trout and mountain whitefish [Prosopium williamsoni]) will influence how Arctic grayling use their habitats. Through this work, I made several contributions towards advancing the scientific exploration of freshwater fish ecology in Arctic watersheds. While freshwater ecosystems compose less than 3% of aquatic habitats worldwide, population decline and extinction in these ecosystems are happening at a disproportionately high rate compared to terrestrial or marine ecosystems (Dudgeon et al. 2006; Reid et al. 2018). The amplification of the biodiversity crisis in freshwater ecosystems is likely to be exacerbated in Arctic watersheds as the Arctic is experiencing warming up to four times faster than the global average (Rantanen et al. 2022). Arctic grayling and their European counterpart, European grayling (Thymallus thymallus), are Holarctic species whose distributions 32 span the high latitudes of North America and Eurasia (Kreuger 1981). Globally, their conservation status is considered Least Concern, though this is changing near the southern peripheries of their distributions where they are being reclassified as Special Concern or at risk in many watersheds (Lashmar and Ptolemy 2001; McPhail 2007; Cahill 2015; Stamford et al. 2017; Alberta Government 2018; Hagen and Gantner 2019). Due to the sensitivity of Arctic grayling to ecosystem changes, they have been proposed as a Boreal indicator species (Cahill 2015). While the work in this dissertation examined the conservation of Arctic grayling in just one southern periphery watershed, it is my hope that these findings can serve as important baseline metrics and techniques for monitoring Arctic grayling across Arctic ecosystems. With the continued increase of remote sensing techniques being applied to questions in freshwater ecosystems, determining the effective scales of research and monitoring of this species will be useful. This work can be a step towards creating a systematic sampling framework that can be applied across the Holarctic range of Arctic grayling to monitor the integrity of freshwater ecosystems in a warming future. Using the advancements I made in this dissertation for both drone sampling in rivers and modeling acoustic telemetry data, scientists working in rivers have further tools available to them for expanding their studies across larger spatial scales. 33 Chapter 2 - Linking fish acoustic telemetry data to spatial covariates in river networks with spatial capture-recapture models Abstract 1. Spatial capture-recapture (SCR) models extend classical capture-recapture models to include spatially-explicit animal locations and environmental covariates. SCR models have been widely employed in terrestrial studies to predict the population size and densities of animals assuming a closed population over a defined area. In this work, I use a relative density formulation to estimate parameters of habitat use from acoustic telemetry data in a branching river network. This is a novel application of SCR to both aquatic systems and acoustic data. 2. Using these models, I found that the relative density of tagged Arctic grayling peaked at water temperatures of 11.3 °C (25th and 75th quantiles 11.1 - 17.1 °C) and showed a positive relationship to the relative density of tagged bull trout during the summer and winter seasons, but not during the spring. Modelled activity centres of tagged Arctic grayling matched well with known seasonal distribution patterns for this population. 3. I found SCR to be advantageous as it: (a) uses the full detection history of the tagged population to link activity centres (average detected locations) to spatially-explicit environmental covariates (analogous to resource selection at the home range scale); (b) can accommodate discontinuous acoustic receiver arrays, which is beneficial in systems where acoustic receiver loss or redeployment is common; and (c) pairs well with existing descriptive approaches to acoustic telemetry analysis which can provide a priori information about the movement patterns of species. 34 4. Due to the nature of acoustic telemetry datasets in which only the tagged population is available for detection, I highlight several specific considerations and assumptions for using this approach: (a) activity centres are fixed within modeling windows by the underlying closed population model, so this method is best applied to discrete, ecologically-relevant timeframes in which fine-scale movements are not the focus; (b) inferences from acoustic telemetry data depict relative (not absolute) densities of only the tagged population; and (c) spatial tagging effort must be defined in the model to ensure that predictions are not merely an artefact of tagging effort across space and time. When applied following these assumptions, this method is broadly useful for aquatic ecologists as it presents a quantitative way to merge large-scale acoustic telemetry datasets with discrete habitat parameters that drive population distributions through time. Introduction Spatial Capture-Recapture (SCR; alternately Spatially-Explicit Capture-Recapture or SECR) is a class of hierarchical statistical models that extends classical capture-recapture models to include both georeferenced animal locations and covariates (Royle et al. 2014). Compared to their predecessors, which do not accommodate the spatial structure of the data, SCR models are used by ecologists to predict the population size and density of animal activity centres within a defined area (Borchers & Efford 2008; Royle et al. 2014). These versatile models can accommodate observation data from any source in which individuals can be uniquely and repeatably identified at discrete locations in space and time, including (but not limited to) camera traps, acoustic surveys, and hair snares (Royle et al. 2014). SCR models have been applied most frequently to data collected from terrestrial organisms moving freely along both the longitudinal and latitudinal axes 35 of the state space (i.e. the collection of pixels that represent all possible locations for activity centres; Royle et al. 2014). Fewer studies have applied SCR to aquatic systems. Those that have done so have typically modeled movements along one functional dimension (upstream and downstream along a single river channel; Raabe et al. 2013; Haydt et al. 2022), with limited extensions into the 2D aquatic environment (Marques et al. 2012; Pirotta et al. 2014). Applications of SCR models to data collected in river networks are even less common (a limited network was used in Murphy et al. 2021) in part because a branching (dendritic) system cannot be easily represented by one or two dimensions. While a state space represented in 1D cannot accommodate branches in the river system, a 2D state space which captures the extent of a watershed will be dominated by pixels of habitat which are non-traversable (i.e. land) for strictly aquatic species. In terrestrial applications, modeling animal movements among dendritic networks has been addressed using ecological distance models that offer more realistic animal movements based on least-cost pathfinding (Sutherland et al. 2015). However, this solution is unnecessarily complex for strictly aquatic species. The ecological distance model allows for terrestrial organisms to limit movement across undesirable habitat but does not completely restrict it, whereas strictly aquatic organisms simply cannot be on land. A more suitable approach to modeling the activity centres of aquatic animals moving within a dendritic state space is offered in Efford’s (2023) secrlinear package for R statistical software (R Core Team 2023). The package, which is a wrapper for the 2D parent package secr (Efford 2023), applies a linear mask approach in which the dendritic network is reduced to a series of 1D lines connected at confluence points within a 2D space, effectively eliminating all land pixels from consideration. From this, a pairwise distance matrix 36 between all wetted pixels is calculated for use in the models (substituting for Euclidean distance; Efford 2023). Acoustic telemetry (hereafter AT) forms the backbone of many large-scale collaborative monitoring programs in both the freshwater and marine environments (e.g., the Ocean Tracking Network, European Tracking Network, the Great Lakes Acoustic Telemetry Observation System, and the Integrated Marine Observing System, among others; Hussey et al. 2015; Hoenner et al. 2018; Hostetter and Royle 2020; Alós et al. 2022). AT datasets are generated by tagging individuals with an acoustic transmitter and releasing them among an array of acoustic receivers. The receiver array then monitors the approximate spatiotemporal locations of each tagged individual by receiving, decoding, and recording signals from transmitters that pass within the detection range of the receivers (Whoriskey et al. 2019). While AT enables monitoring of aquatic organisms at spatial scales that would not be feasible using other sampling methods, it cannot detect new, untagged individuals within the vicinity of receivers (as opposed, for example, to camera traps used to study individuals with uniquely identifiable markings, e.g. Dorazio and Karanth 2017). For this reason, SCR models fit to AT data alone cannot estimate population size, though they can be suitable for modeling the relative density of a tagged population in relation to spatial covariates. In this paper, I develop SCR models that are suitable for use with AT data in a dendritic river network. I use these models to quantify the relative density of tagged individuals in relation to spatial covariates (which, to aquatic ecologists, is conceptually analogous to resource use, though typical selection analyses don’t consider relative density; sensu Johnson 1980). Through a case study, I demonstrate how quantifying these relationships with spatial covariates can be used to explore ecological hypotheses. Specifically, using the packages secr and secrlinear (Efford 2023), I assess whether the seasonal distributions of Arctic grayling (Thymallus arcticus, Pallus 37 1776) in the Parsnip River watershed of north-central British Columbia are influenced by two spatial covariates: watershed-scale temperatures and the relative density of sympatric bull trout (Salvelinus confluentus, Sukley 1859). As high temperatures likely limit the distributions of summer rearing Arctic grayling (Stamford et al. 2017), I expect to find a negative response between high temperatures and the activity centres of tagged Arctic grayling. If the trophic relationship between Arctic grayling and bull trout is antagonistic, it is likely that Arctic grayling foraging activity would be reduced in cases where bull trout have high relative densities, and I expect a negative association between the two. I conducted my analysis across three ecologically-relevant periods in which adult Arctic grayling use spatially discrete habitats for overwintering, spring spawning, and summer rearing (Blackman et al. 2002; Stamford et al. 2017). Methods A primer for SCR models and their application to AT in dendritic networks The main data collected from AT are time-stamped detections from n individuals at J receivers. Specifically, SCR models ‫ݕ‬௜௝௞ , which denotes whether the tagged individual ݅ was detected by acoustic receiver ݆ during occasion ݇. For this analysis, occasions represented aggregated detections over each day. The first component of an SCR model is the observation function, which represents the imperfect detection of tagged individuals by the receiver array. For AT receivers, the most appropriate observation model (alternately, ‘detection function’) offered is the ‘proximity’ type, wherein each receiver can observe multiple unique tags in a day and individuals are not precluded from being detected at other locations after their first detection. The proximity detector type differs from the more literal capture-recapture methods which do not allow for multiple spatial encounters until an animal is released (Efford et al. 2009; Royle et al. 2014). In the observation 38 model, the probability of detecting an individual decays as a function of the distance between the location of an individual’s latent activity centre ࢙௜ , and that of the receiver recording a detection, ࢞௜ . Specifically, I used a half-normal distance function: మ ିฮ࢞ೕ ି࢙೔ ฮ ݃൫࢙௜ , ࢞௝ |ࣂ൯ = ߙ exp ቆ ଶఙమ ቇ, (Eq. 2.1) where ݃൫࢙௜ , ࢞௝ |ࣂ൯ is the probability that receiver j detects an individual i during an occasion given the vector of detection parameters ࣂ; ฮ࢞௝ − ࢙௜ ฮ is the distance between the latent activity centre i and the acoustic receiver j; ߙ is the probability of detection when an individual’s centre of activity is located exactly at the location of the receiver; and ߪ is a scale parameter whose magnitude determines how detection probability decays as a function of distance between activity centres and the receiver (the larger ߪ, the slower the decline in detection probability with distance). To use this SCR observation model with data collected by acoustic receivers in a dendritic river network, one simply replaces the Euclidean distance term ฮ࢞௝ − ࢙௜ ฮ with network distance, or the pairwise distance matrix between pixels across the linear state space (Efford et al. 2009); although, more sophisticated distance terms have been developed (e.g. Sutherland et al. 2015). In practice, I am modelling an aggregate of the binary encounters of an individual across all receivers, known as the detection history, ࢟௜ = (‫ݕ‬௜ଵ , … , ‫ݕ‬௜௝, … , ‫ݕ‬௜௃ ), where ‫ݕ‬௜௝ , represents the number of occasions when individual i was detected at receiver j. If I assume that the detection function (Eq. 2.1) remains constant across all K occasions, the probability of observing the detection history of an individual, ࢟௜ , given its activity centre was at ࢙ is defined as: ௃ Pr(࢟௜ |ࣂ, ࢙) = ∏௝ୀଵ Binomial ቀ‫ܭ‬, ݃൫࢙, ࢞௝ ൯ቁ (Eq. 2.2) 39 The second component of the SCR model represents the relative densities of the latent activity centres of tagged fish as a function of some set of ecological, environmental or sampling covariates. Specifically, it models the location of the activity centres as an inhomogeneous Poisson point process, where the density of activity centres at location s, ߤ(࢙), is a function of V spatial covariates. As in a generalized linear model, I can ensure that the density values remain positive by using a natural log link function: ln൫ߤ(࢙|ࣘ)൯ = ߚ଴ + ∑௏௩ୀଵ ߚ௩ ܿ௩ (࢙) + ߚா ܿா (‫)ܛ‬, (Eq. 2.3) where ߚ଴ is the intercept, and ߚ௩ is the coefficient describing how the density of activity centres changes with values of covariate v at location s, denoted by ܿ௩ (࢙). In addition to the spatial covariates of ecological interest, when using AT datasets it is important to include a covariate that represents tagging effort, ܿா (‫)ܛ‬, and its coefficient, ߚா . Including this spatial covariate for tagging effort is necessary to ensure that the predicted relative density is not simply reflecting tagging effort (defined in my case study as the cumulative kernel density estimate of tags applied at each pixel for each season). Many of the relative density parameters ࣘ = (ߚ଴ , ߚଵ , . . . , ߚ௏, ߚா ) represent the relationships with a spatial covariate and estimating them is an important inferential target of the method. Incorporating this point process model for AT data requires modifications. Unlike most SCR models in which the objective is to estimate the size and absolute density of the population, here I can only detect a fixed number of tagged animals n. Thus, I simplify the maximum likelihood function described in Borchers and Efford (2008) to one that maximizes the likelihood of the detection and density parameters conditional on n to predict the relative density of tagged animals in space. The relative density model is defined as (Efford 2023): 40 ݊ ℒ(ࣘ, ࣂ|࣓) = ቀ݊ … ݊ ቁ ∏௡௜ୀଵ ∫ோమ Pr(࢟௜ |ࣂ, ࢙) ߤ(࢙|ࣘ)࢙݀, ଵ ௖ (Eq. 2.4) where ߤ(࢙|ࣘ) is the distribution of tagged animals in the state space as a function of spatial covariates; Pr(࢟௜ |ࣂ, ࢙) is the probability of observing the detection history ߱ for individual ݅ given its activity centre was at ࢙; (Efford 2023; Borchers and Efford 2008). For a more detailed description of the likelihood functions used by secr and secrlinear, please refer to Efford et al. (2009) and Borchers and Efford (2008). Study Area The data used in my case study were collected in the Parsnip River watershed (approximately 5,000 km2) in north-central British Columbia, Canada (Figure 2.1). The area is located on the overlapping traditional territories of Treaty 8 signatory Nations (Prophet River, Saulteau, and West Moberly First Nations) and the traditional and unceded territory of the McLeod Lake Tse’Khene. The Parsnip River is a moderate- to low-gradient system which flows northwest along the Rocky Mountain Trench. Over its 240 km course, it is fed by many low-order tributaries draining the Hart Ranges as well as four higher-order rivers which support Arctic grayling and were selected to be included in this study: the Anzac, Table, Hominka, and Missinka Rivers. Together, the Parsnip River and these main tributaries contain most of the Arctic grayling in the Parsnip River watershed (Stamford et al. 2017). Acoustic tagging Adult Arctic grayling and bull trout were captured by angling with either dry flies or spin casting based on current angling conditions and hatches. Tagging began after spring freshet and continued 41 until ice-up at various sites in the Parsnip River and its tributaries. Captured Arctic grayling and bull trout > 230 g were measured for fork length (mm) and mass (g) before being surgically tagged with acoustic transmitters (Innovasea [formerly Vemco] V9 tags, 90-150 sec transmission interval, Bedford, NS; full tagging procedures are available in Supplemental information). Fish were immobilized with electroimmobilization gloves (Smith-Root, Vancouver, USA) for surgeries and only individuals large enough to accommodate tag burden (> 230 g) were used in this study. All individuals were captured as authorized under permits PG18-356580, PG19-523435, PG20606121, and PG21-622265 issued by the BC Ministry of Forests, Lands, Natural Resources Operations and Rural Development (currently BC Ministry of Forests). Fish handling and tagging were conducted in accordance with the University of Northern British Columbia’s Animal Care and Use Committee (ACUC protocol 2018-06). To distribute capture effort across the watershed, I followed a structured approach with approximately equal tagging effort in the upper, middle, and lower reaches of each tributary (with some additional effort applied on the mainstem of the Parsnip River) during the summer when adult Arctic grayling position themselves in tributary feeding areas (McPhail 2007, Hughes 1999). Given the heterogeneous distribution of fish across the watershed through the summer period, tag distribution was ultimately unequal between different sections of the watershed (Table 2.1). Data Collection Innovasea VR2W- and VR2Tx-69 kHz acoustic receivers were deployed over an extent of 366 river kms. The selection of deployment sites was based on site accessibility, substrate conditions (coarse enough to not lose receivers to silt compaction while fine enough to minimize acoustic interference from water flowing over large substrates), distance from public access points (accessible but not likely to be seen by anglers), qualitative habitat suitability for Arctic grayling 42 (areas of flow refuge and forage opportunities during summer feeding), and amenability to receiver redundancy within 2 km. Receivers were deployed following a clustered strategy (2 – 4 receivers spaced 0.5 – 2 km apart) while distributing the receiver clusters widely in the study rivers to optimize data collection for use with SCR models (Royle et al. 2014). The clustered deployment also offered redundance against receiver loss or damage due to freshet, scouring debris, shifting substrates, or vandalism. The receivers were moored hydrophone-up at each site using a cement block, steel cable, and a cable anchoring system (Duckbill Earth Anchors, Fort Mill, SC). Deployment of the acoustic receiver array began in late July 2018, with 54 receivers deployed across the watershed by the end of the open-water season. Deployments and lost receiver replacements continued after freshet in 2019 (27 lost; 41 deployed; 68 total), 2020 (18 lost; 2 deployed; 52 total), and 2021 (12 lost; 10 deployed; 50 total) (Table 2.2). Deployed receivers were downloaded and maintained once per year during autumn base flow conditions. Equipment retrieval began in late summer of 2021 following approximate downriver migration timelines and the monitoring window officially concluded on 30 October 2021. Fitting SCR models to detection data SCR models can accommodate spatial covariates on the condition that they are available across the state space used to model activity centres (with additional flexibility to integrate receiver-level covariates and individual-level covariates). The models in my case study used a dendritic statespace of 366 pixels at a resolution of one river kilometer (Rkm) and covariates were defined to the same resolution. Two spatial covariates were included in my SCR model of Arctic grayling: water temperature and relative density of bull trout. Temperature for each state-space pixel was obtained using predictions from a spatial stream network model parameterized with data collected from 221 temperature loggers installed throughout the watershed over the same spatiotemporal extent as the 43 acoustic receiver array (O’Connor 2023). The temperature covariate was used in the summer rearing period models (when high temperatures have the potential to be a limiting factor to Arctic grayling distributions) and represents the mean pixel predictions of daily temperatures (Temp) from July 1 – September 15. The second covariate was a relative density surface of tagged bull trout activity centres throughout the watershed (BT), as characterized by an intercept-only SCR model for bull trout detection data collected across each sampled window (Supplemental Information SI 2. 5). Since the Temp and BT covariates had different magnitudes, the covariates were standardized for analysis (Schielzeth 2010). Due to the sensitivity of presenting spatial information on the distributions of blue-listed bull trout in the province of British Columbia, spatial plots of the BT covariate are not presented here but may be made available upon reasonable request. Spatial tagging effort is also included as a covariate in the model. Three sets of Arctic grayling SCR multisession models were fit for the spring, summer, and winter detection windows. The summer models (Jul. 1 to Sep. 15) were selected from six candidate models, and the spring (May 1 to Jun. 1) and winter (Nov. 1 to Mar. 31) modeling windows were each selected from two candidate models (Table 2.3). The dates used to bound seasonal detection windows were defined as core periods during each seasonal life history stage (overwintering, spring spawning and summer feeding). Shoulder seasons during which fish move and transition between life history stages were omitted from the analyses. For example, the spring spawning migration was defined as the month of May, while in practice it can extend into June and varies with ice-out timing (Blackman 2002). Conversely, early summer rearing can begin as early as June and the transition between the spring spawning and summer rearing periods can be relatively fluid from a data bounding perspective. No models were fit to 2018 data as both the acoustic receiver 44 array and tagging program were in their first-year deployment phases and detections were considered preliminary as the array coverage was incomplete and the sample size was low. For the summer feeding period, candidate models included the no effects model (i.e. no covariates, or the intercept-only model of activity centre density), models including standardized temperature as both linear and polynomial predictors (Temp and Temp + Temp2; the polynomial term is included as there are likely upper and lower limits to thermal habitat use in this watershed), a model with standardized relative tagged bull trout activity centre density (BT), and models with all combinations (Temp + BT and Temp + Temp2 + BT). For the spring and winter windows during which temperature predictions were not available, the no effects model was compared against the model including only the tagged bull trout relative density covariate BT. Models for each season were fit as multi-session models pooling seasonal data from 2019, 2020, and 2021 with annually defined spatial covariates varying by session, but assuming common effects of Temp, Temp + Temp2 and BT across sessions. Model selection was conducted using Akaike’s Information Criterion with small sample correction (AICc; Burnham and Anderson 2002; Table 2.3). Model predictions were averaged based on the 95% confidence set for the top models selected by AICc (i.e. by averaging predictions of models with cumulative AICc weights that added up to at least 0.95; Burnham & Anderson 2002). Temperature predictions were made with the standardized BT covariate set at zero, which was the equivalent of the mean of the untransformed variable. The codes for fitting and selecting models as well as extracting and visualizing activity centres are available in Appendix A. 45 Results The acoustic receiver array recorded over three million acoustic detections which were aggregated into 5771 unique daily Arctic grayling locations. Model selection indicated D ~ Temp + Temp2 + BT as the best model for the summer multisession model of relative density D, with secondary support for model D ~ Temp + Temp2 (Table 2.3; Supplemental Information SI 2.1). The model D ~ BT was the best model for the winter season, and the no effects model was selected as the best model for the spring season with some uncertain effect of the D ~ BT model indicated by AICc weight (Table 2.3; Supplemental Information SI 2.1). Temperature predictions for tagged Arctic grayling activity centres peaked at 11.3 °C with approximately 1.7 tagged individuals per Rkm (25th and 75th quantiles at 11.1 and 17.1 °C; Figure 2.2a). The density of tagged Arctic grayling increased with the density of tagged bull trout in the summer and winter, with the strongest relationship occurring in the winter (Figure 2.2b). During the summer feeding window, tagged Arctic grayling were widely distributed throughout the Parsnip River and the four major tributaries included in this study (Figure 2.3). High relative densities were estimated to reach the uppermost reaches of all study tributaries, and notable densities occurred in the Parsnip River mainstem near the Table River confluence. The model also revealed more apparent spatial separation between activity centres in the summer than during other seasons (Figure 2.3, Figure 2.4, Figure 2.5). During spring spawning, Arctic grayling relative densities were distributed widely in the Parsnip River mainstem and in the lower reaches of all four study tributaries with a notably high relative density in the Parsnip River just north of the Table River confluence. Spring distributions in the Table River were also centred in middle and upper reaches, higher than in the Anzac, Hominka, 46 and Missinka Rivers (Figure 2.4). During winter, activity centres for Arctic grayling activity were predominantly confined to the Parsnip River mainstem with notable concentrations in the lower reaches of the Anzac River and Table Rivers (Figure 2.5). Most overwinter activity centres were proximate to the tributary confluences of the Anzac and Table Rivers with the Parsnip River mainstem. Models were predicted using the summer temperature distributions in a cool year (2019) and a warm year (2021) to create a watershed-scale estimate of Arctic grayling densities as a function of temperature (Figure 2.6). During the cool year, activity centres were widely distributed across the tributaries. These ranges contracted during the warm year into a patchier distribution. The Parsnip River mainstem was predicted as low density in both years, though some habitat use that was predicted to be used in the cool year contracted in the warm year. Discussion This study demonstrated the usefulness of SCR as a quantitative tool for analysing AT data in dendritic river networks. An attractive feature of AT data for aquatic ecologists is its ability to describe where animals distribute themselves through space and time. A large proportion of AT studies to date have focused on descriptive approaches, and fewer have modeled the drivers of animal density and distribution (Whoriskey et al. 2019; Brownscombe et al. 2022). Formalizing an approach to AT analysis that quantifies relationships between spatial covariates and tagged animal density could be a step towards unifying the vast amount of data produced by AT studies and the quantitative objectives which often drive resource management (Crossin et al. 2017). The SCR models that I developed were useful for describing how a small set of spatial covariates influenced the watershed-scale distribution of tagged Arctic grayling. The response of tagged fish 47 to temperature showed a peak of 11.3 °C (25th and 75th quantiles 11.1 – 17.1 °C), which matched well with the results of other thermal habitat studies. Habitat suitability indices developed in other systems for Arctic grayling suggested optimal temperature ranges of 5 – 14 °C (all life stages; Larocque et al. 2014) and 6 – 16 °C (adults; Hubert et al. 1985; Stewart et al. 2007). The sharp decline in the response at temperatures above 16 °C broadly agrees with reduced survival rates above 18 °C identified in Arctic grayling populations in the southern periphery population of Big Hole, Montana (Carillo-Longoria 2023). A habitat occupancy study in the Parsnip River watershed found a high probability of occupancy by Arctic grayling in this system at 10.9 °C (8.7 – 14.2 °C; O’Connor 2023). A series of thermal preference experiments conducted on adult Anzac River Arctic grayling found a thermal preference peak at 11.6 °C (25th to 75th quantile range of 10.1 – 13.0 °C; Bottoms et al. unpublished data). These findings also reflect a selection for cooler temperatures found in a study of habitat use by Arctic grayling in an in-situ study conducted during the same year as the thermal preference work (Bottoms et al. unpublished data). Maps of the spatial distribution of activity centres and density of tagged fish showed distinct seasonal patterns, which matched what has been broadly described in this watershed by other studies (Blackman 2002; Stamford et al. 2017; Hagen and Stamford 2017; Supplemental Information SI 2.3). Model outputs indicated persistent use of overwinter habitat in the lower Anzac and Table Rivers (in addition to the well-established winter habitats in the Parsnip River mainstem) which was first documented in a radio tagging study by Blackman (2002). I note that the assumption of the closed population model implies that activity centre locations are fixed in each season (Efford 2023). This assumption is reasonable for my study, as seasonal distributions are well-understood and Arctic grayling are thought to remain close to their territorial feeding sites during each modeling window. 48 While I conceived models with the BT covariate (Supplemental Information SI 2.2) to investigate Arctic grayling distribution in response to the relative density of as a potential predator, the precise trophic relationship between the two sympatric species remains unknown (Stamford et al. 2017). A negative relationship (avoidance) with the bull trout covariate could have been evidence to support a predatory relationship, which in turn may have warranted considering them as a limiting factor in the conservation plan for Arctic grayling in the Parsnip River watershed (Stamford et al. 2017). Instead, I found a positive relationship between Arctic grayling and bull trout during the summer and winter seasons. A positive relationship could result from a predatory behavior, where bull trout actively seek out areas of Arctic grayling occurrence, or it could be a function of shared habitat needs in which co-occurrence is made possible through adaptive behaviors (e.g. foraging at different niches within the same habitats; sensu Nakano 1998). If the relationship is indeed predatory, it is likely that Arctic grayling foraging activity is reduced in cases where bull trout have high relative densities. A stable isotope study conducted on this population found that bull trout feed on prey at the same trophic level as Arctic grayling (Clevenger et al. unpublished data), and while the size of adult Arctic grayling is often below the gape limitation threshold of bull trout (Stamford et al. 2017), evidence in other systems suggest that Arctic grayling may use their prominent dorsal fin to appear larger and deter predation (Stewart et al. 2007). There was no support for the bull trout covariate in the top spring model, though the model with the bull trout covariate had some support based on AICc weights. This could likely result in a mismatch in migration timings where adfluvial bull trout enter the tributaries from the Williston Reservoir later in the summer months and overlap was not likely to occur in many cases (Hagen and Weber 2019). The detection data and corresponding density surface predictions for the BT covariate were relatively sparse compared to the Arctic grayling dataset, as only 42% of the tags in the study were 49 allocated to bull trout and a further portion of these individuals outmigrate from the study system to overwinter in the Willison Reservoir (Langston and Cubberley 2008; Stamford et al. 2017; Supplemental Information SI 2.4). While this study demonstrates that examining interspecific relationships can be achieved with these methods, it may be that covariates based on co-occurring species may derive a better signal relevant to behavioural studies when focused on species with well-defined trophic niches and predator responses. When applying SCR to AT data, only tagged individuals can be detected, and it is important to consider whether the sampling design is appropriate and to verify that what is being modeled represents the distribution of tagged animals rather than a distribution of tagging effort. While I attempted to apply a structured tagging approach across streams, heterogeneous distributions of site conditions and both study species over time led to some unequal distributions of tags within the array. I compensated for this by including tagging effort as an external spatial covariate when evaluating the likelihood (Efford 2023). In my case study, I only had data on when and where tags were deployed. A better metric to be considered for future AT studies would include hours spent tagging at sites with both high and low catch rates to compute a metric of catch per unit effort. I further compensated for my variable coverage by modeling my covariate predictions as multisession models with seasonally-specified array deployment files indicating whether each receiver was online or offline on each occasion. Receiver loss was prominent in the 2019 and 2020 freshets and sustained high-water years (particularly in the Table and Anzac Rivers; Supplemental information SI 2.6), and while suitable sites were replaced, AT detections were reduced in the spring and summer in these tributaries. As a function of applying AT in a riverine system with highly variable discharge, my detection efficiency was also low at measured sites during high flows (a consideration of AT use in rivers in 50 general). The data included in my model were filtered for false and out-of-bounds detections, though skip movements (cases where a tag was detected at two receivers but missed by a receiver located between the two) were permitted to retain valid detections which would otherwise be filtered out because of variable detection efficiency. I compared the model outputs in my case study to models fit to the same dataset filtered for skip movements with the R package actel (Flávio & Baktoft 2021) and both model performance and predictions were similar. I chose to use the dataset which permitted skip movements as it reduced the number of activity centres being predicted in seasonally unlikely places in cases where a tagged individual rapidly migrated up or downriver near the peripheries of the monitoring periods but was missed by the receiver array. I recognize that many formulations of SCR models are designed for population estimation and I reiterate how AT data are not suited for this objective. Indeed, this formulation is overly complex for my needs and my use of the relative density function in secr simplifies the underlying model. I see further potential for the development of specific models for use with AT data. For example, Raabe (2013) presented a Bayesian formulation of an SCR model that could be a suitable candidate to be modified for use in river networks following Efford’s (2023) approach to replacing the Euclidean distance term d with network distance. I also see high potential for further development of this class of models to be useful in analyzing AT data with other spatial covariates interpolated over various scales (e.g., data collected via remote sensing by satellites or drones). Assumptions of this approach are as follows: (i) this formulation works using a closed population model and assumes static activity centres over each modeling window; (ii) insights derived from this approach apply only to the tagged population and careful consideration of a species’ ecology should be used when extrapolating these findings across unsampled units; (iii) tagging effort should be distributed evenly over the sampling area and a spatial layer of tagging effort should be 51 applied to compensate for any residual unevenness in tagging; and (iv) tagged individuals are assumed to survive over the duration of the tagging window. While I apply this approach to a dendritic river network, the underlying assumptions that apply to AT data in river networks would also apply to AT data collected in lakes, estuaries and oceans, though careful consideration should be given to how 2D state spaces differ from my case study. The assumptions of a closed population and static activity centres over each modeled window would also apply to open-water species provided that a researcher can define ecologically-relevant periods (e.g., seasons) where the individuals have settled in static home ranges. Computation times, which in my dendritic state space were relatively efficient (~5 minutes for the most complex models on a mid-range machine), would increase in proportion to the extent of the area represented by a 2D state space, which while not prohibitive, warrants consideration among users when specifying the spatiotemporal grain and extent of their covariate pixels. In conclusion, this work demonstrated the use of a powerful quantitative approach for analyzing the vast amount of AT data produced by fisheries ecologists. I found SCR modeling of AT data to be useful and suggest considerable potential for this class of modeling in helping bridge the gap between traditional telemetry analyses and the quantitative objectives that often drive management and conservation efforts (Crossin et al. 2017). I found that combining the statistical model with existing AT visualization techniques (e.g. animated tracks produced with the R package RSP; Niella et al. 2020) produced the most comprehensive understanding of Arctic grayling distribution within the available spatiotemporal dimensions of the dataset. I submit that SCR models have strong potential for further development with a specific focus on AT data. While movements between seasons are implied, they are not explicitly modeled in my case study, and further development of movement models which use the scale parameter ߪ with respect to spatially 52 varying environmental gradients may prove useful (Royle et al. 2014). Through modifications of the detection function, there may be opportunities to integrate independent data that reflect detection efficiency. At the time of writing, there is promising ongoing development for richer analyses in the pursuit of applying SCR models to AT data (Hostetter and Royle 2020), including generalized telemetry accommodations in the secr package for non-linear state-spaces (Efford 2024, in prep). 53 Totals 2021 2020 2019 2018 Year Reach Lower Middle Upper Lower Middle Upper Lower Middle Upper Lower Middle Upper 8 1 7 12 2 2 8 Parsnip GR BT 6 52 10 7 9 7 1 4 2 Anzac GR BT 20 50 1 11 5 8 12 3 10 24 1 14 3 3 1 2 Table GR BT 24 3 6 15 2 14 3 5 4 Hominka GR BT 13 3 10 5 2 1 2 Missinka GR BT Totals GR BT N 20 0 20 10 9 19 19 4 23 0 5 5 33 8 41 19 9 28 0 2 2 18 4 22 21 18 39 0 0 0 1 1 2 6 2 8 147 62 209 54 Table 2.1. Annual counts of Arctic grayling Thymallus arcticus (GR) and bull trout Salvelinus confluentus (BT) acoustically tagged in the Parsnip River and its four major tributaries, British Columbia, Canada. The species codes displayed correspond to the codes used by the Ministry of Environment for British Columbia. Chapter 2 Tables Parsnip 28 15 1 1 Year 2018 2019 2020 2021 5 1 3 11 Anzac 3 0 2 13 Table 1 0 11 1 Hominka 0 0 10 1 Missinka 50 (+10 / -12) 52 (+2 / -18) 68 (+41 / -27) 54 Total 55 Table 2.2. Annual deployments and active coverage extent of the acoustic receiver array in the Parsnip River and its four major tributaries, British Columbia, Canada. High receiver loss occurred during a high discharge freshet in 2020. See Figure 2.1 for a visual description of these data. Values in parenthesis represent the receiver deployments (+) and losses (-) during each year. 1340.24 1346.43 1340.95 1346.84 4 3 D ~ BT D ~ 1 (No effects) -666.12 -670.21 AICc Model AICc 3205.08 3207.00 3213.58 3214.27 3223.09 3228.46 Winter models k logLik AIC k 6 5 5 4 4 3 AICc 1107.20 1109.45 Model D ~ 1 (No effects) D ~ BT Model D ~ Temp + Temp2 + BT D ~ Temp + Temp2 D ~ Temp + BT D ~ Temp D ~ BT D ~ 1 (No effects) Summer models logLik AIC -1595.81 3203.63 -1597.99 3205.98 -1601.28 3212.56 -1602.80 3213.60 -1607.21 3222.42 -1611.03 3228.07 Spring models k logLik AIC 3 -550.38 1106.76 4 -550.36 1108.71 0.00 5.89 dAICc dAICc 0.00 2.26 dAICc 0.00 1.92 8.50 9.19 18.01 23.38 0.95 0.05 AICcwt AICcwt 0.76 0.24 AICcw 0.71 0.27 0.01 0.01 0.00 0.00 56 Table 2.3. AIC statistics for candidate Arctic grayling multisession SCR models for each season. Presented are the log likelihood, the standard AIC score and the AICc score for small sample sizes, the difference in AICc between each model from the top model (ΔAICc), and the respective model weights (AICcw). Model terms are the relative density response (D), temperature and squared temperature, and the bull trout covariate BT. Chapter 2 Figures 57 58 Figure 2.1. The Parsnip River watershed in north-central British Columbia, Canada, was discretized into a state-space mask of 366 pixels at a resolution of 1 Rkm/pixel. The state-space mask includes the Parsnip River mainstem and four major tributaries included in this study: the Anzac, Table, Hominka, and Missinka Rivers. The watershed flows northwest through the Rocky Mountain Trench into the Williston Reservoir. Acoustic receiver deployments are indicated by year; overlapping points indicate replacements of compromised or lost acoustic receivers. 59 Figure 2.2. Model-averaged predictions of the relative density of tagged Arctic grayling activity centres with relation to the temperature covariate across all three summer seasons (panel A) and the bull trout covariate BT in the spring and summer (panel B). Relative density of tagged Arctic grayling peaked at 11.3 °C (25th and 75th quantiles 11.1 – 17.1 °C). The shaded region indicates the ± 1 standard error of model predictions. Tagged Arctic grayling relative density response to tagged bull trout relative density was positive in both the winter and summer, with a stronger response in the winter. Areas where the distributions of standard error width is narrow have few data points of activity centre density estimates at those temperatures in 2019 and 2020. 60 Figure 2.3. Annual summer distributions of tagged Arctic grayling activity centres estimated by the SCR model across the Parsnip River watershed. Activity centres are given by white markers over the relative tagged population density across the state-space in blue. Coloured circles indicate the receiver array and which acoustic receivers were active over the detection window (Jul. 1 – Sep. 15 each year). Summer activity centre distributions were wide ranging within the Parsnip River and its four main tributaries with mainstem Parsnip River activity centres being closely associated with the Table River confluence. 61 Figure 2.4. Annual spring distributions of tagged Arctic grayling activity centres estimated by the SCR model across the Parsnip watershed. Activity centres are given by white markers over the relative tagged population density across the state-space in blue. Coloured circles indicate the receiver array and which acoustic receivers were active over the detection window (May 1 – Jun. 1 each year). Spring activity centre distributions occurred largely in the Parsnip River mainstem and the lower to middle reaches of the tributaries. 62 Figure 2.5. Annual winter distributions of tagged Arctic grayling activity centres estimated by the SCR model across the Parsnip watershed. Activity centres are given by white markers over the relative tagged population density across the state-space in blue. Coloured circles indicate the receiver array and which acoustic receivers were active over the detection window (Nov. 1 – Mar. 31 each year). Winter activity centre distributions were concentrated to the Parsnip River mainstem and the lower reaches of the Anzac River. One activity centre appearing in the upper Missinka River in 2021 is suspected to be a case of mortality as the tag ID does not reappear in the subsequent spring or summer models for 2021. 63 Figure 2.6. Model-averaged spatially-explicit predictions of the density of Arctic grayling activity centres in the Parsnip River watershed based on the mean temperature at each location in the state space in the 2019 (a warm year) and 2021 (a cool year) summer seasons. Distribution is predicted to be higher in the tributaries and upper Parsnip River where summer temperatures are cooler than the Parsnip River mainstem. Predicted densities show wide Arctic grayling distributions during the cool year and become patchier and more discrete during the warm year. Chapter 3 - The distribution of Arctic grayling (Thymallus arcticus) is related to temperature at the reach scale, but not at the river scale Abstract 1. As a result of the continuous nature of river habitats, the drivers of the abundance and distribution of riverine fishes are often multiscale. Studies within a single scale can miss important environmental predictors of these drivers and often leave conservation managers without information critical to the recovery of a species. Studies that apply a hierarchical, nested study design can fill these gaps and identify environmental predictors which apply across scales. 2. To demonstrate the importance of studying habitat use at different scales, I investigated the use of physical and thermal habitats by Arctic grayling at two spatial scales within the Anzac River of north-central British Columbia, Canada. As an index of their abundance, I counted Arctic grayling using snorkel surveys during their summer feeding period at both the reach (100 m – 2 km) and river (1 – 52 km) scales. I used a visual/thermal drone to create metrics of the physical and thermal habitats available to the surveyed fish. 3. I used generalized linear (mixed) models to understand how Arctic grayling abundance indices are affected by environmental and biological covariates. Pool habitat was an important predictor of Arctic grayling abundance at both the reach and river scales. Reachscale models revealed that abundance was greater in larger pools, and river-scale models revealed that abundance was positively related to both the mean ranked size and density of pools within each 1-km sampling unit. In contrast, temperature was only an important predictor at the reach scale. The relationship between the Arctic grayling abundance index and cool temperatures at the reach scale matched the one found at the watershed scale in 64 Chapter 2. Unexpectedly, upstream distance superseded temperature as an environmental predictor of abundance index at the river scale. 4. Conservation planning for rivers is often focused on preserving coldwater-producing habitats (i.e. headwaters). My findings support the results of other multiscale investigations in rivers and highlight the important effects of non-headwater habitats which, while warmer, support varied instream habitats and may be important to conservation actions through land procurements and protections. Introduction Determining the appropriate scale to conduct studies linking fluvial fish ecology to conservation has been a recurring theme in the riverine literature (Fausch et al. 2002; Torgersen et al. 2022). Rivers are relatively narrow corridors of habitat that cut, meander, or braid across the larger landscape (i.e. riverscapes; Fausch et al. 2002). The processes that drive the abundance and distribution of organisms within rivers are often multiscale – that is, they are dependent on both the reach-scale heterogeneity of the rivers themselves (Hughes 1998; Sears et al. 2019) and the context of the large-scale landscapes through which they flow (Morash et al. 2021). Studies at the river scale (10s – 100s of km) can be useful for broad watershed management, though important, and often overriding, effects of discrete habitat features (e.g. thermal refugia associated with groundwater upwelling) can be lost at this coarse resolution (Fausch et al. 2002; Dzara et al. 2019). Studies at the reach scale (100s – 1,000s of m) can both capture the effects of these features and be useful for studying fish behaviour or physiology, but they often lack the context of, and thus applicability to, the river at large. Rivers themselves also demonstrate processes at multiple temporal scales, exhibiting variability in the short-term (e.g. daily thermal cycles), intermediate65 term (e.g. seasonal, hydrological and thermal cycles), and, in some cases, the long-term (e.g. annual and interannual shifts in the course of the river itself) (Caissie 2006; Stanford et al. 2005; Hauer et al. 2016). Collecting comprehensive fisheries data in riverine systems can be logistically challenging for a myriad of reasons related to topography, access, hydroperiod, migration timing, and spatiotemporal scales of interest (Fausch et al. 2002). Consequently, most fisheries studies have been conducted at intermediate scales over which study objectives can be effectively explored. Inferences from these scales are extrapolated across unsampled units to estimate ecological patterns across larger scales (Torgersen et al. 2022). Rivers are spatially continuous habitats that demonstrate longitudinal connectivity in both physical structures and physicochemical properties along a gradient from headwaters to mouth (Vannote 1980; Caissie 2006). As such, unit sampling approaches can be limited in capturing the full variability of habitat features among all sites within a river (Schlosser 1991; Poole 2002; Torgersen et al. 2022). To overcome these limitations, Fausch et al. (2002) proposed a riverscape sampling design that applied a hierarchical, nested approach to collecting fisheries data, thereby contextualizing reach-scale inferences into the large-scale characteristics of riverine systems. Studies that have employed riverscape sampling approaches have revealed the importance of multiscale effects on ecological patterns and processes in stream fishes. In the Arkansas River basin, variation in the abundance of Arkansas darter (Etheostoma cragini) was explained by canopy cover and channel depth at the reach and segment scales, but at the basin scale the best environmental predictors were stream width and depth (Wellenmeyer et al. 2019). White-spotted char (Salvelinus leucomaenis) and masu salmon (Oncorhynchus masou) showed multiscale movement patterns where body condition was found to drive emigration behaviours across local 66 and tributary scales in a Japanese stream system (Kanno et al. 2021). Local-scale distribution of pool habitats and thermal refugia were found to be more important than river-scale distributions of temperature for Chinook salmon (Oncorhynchus tshawytscha) in an Oregon watershed (Torgersen et al. 1999). While applications of riverscape sampling have increased since the publication of Fausch et al. (2002), few studies have explicitly addressed how inferences derived from this approach can improve conservation management (Torgersen et al. 2022). Indeed, problems in conservation management are also multiscale and the mismatch between scales used by scientists and conservation planners has been identified as a barrier to implementing effective conservation policy (Guerrero et al. 2013; Rose et al. 2018). The conservation of Arctic grayling in the Williston Reservoir watershed of north-central B.C is an example of the mismatch between the scale of research and conservation. The seasonal habitat use of Arctic grayling (Thymallus arcticus) occurs over discrete (reach) scales, and their migrations between these habitats occurs over river and watershed scales (Blackman et al. 2002). It is at the river scale where research has been primarily focused (e.g. Blackman et al. 2002; Hagen and Gantner 2019). However, choosing specific sites for conservation action through boots-on-the-ground interventions is the current focus of conservation managers in this system (e.g. Stamford et al. 2017), and this requires contextualizing fine-scale research across the greater watershed to highlight discrete areas where specific actions can have the most benefit. In the Williston Reservoir watershed, freshwater habitats have been extensively altered following the completion of the W.A.C. Bennett Dam in 1967. The conservation of riverine species in this watershed requires knowledge of both the drivers of their habitat use (e.g. temperature, instream habitats) as well as how these drivers are distributed across the landscape. The effective use of 67 conservation strategies in this watershed will further require knowledge of how these drivers differ across spatial and temporal scales. Populations of Arctic grayling have been particularly impacted by habitat loss in this watershed as the flooding of the Williston Reservoir converted extensive areas of riverine habitats into slower, deeper, more lake-like habitats (Lashmar and Ptolemy 2002). The effects of habitat loss were compounded by overfishing and industrial development leading to the imposition of a harvest moratorium in 1995 and the ongoing development of a conservation plan for the species (Blackman 2001; Stamford et al. 2017). Through the development of this conservation plan, summer water temperatures were identified as a probable limiting factor preventing the recovery of the species (Stamford et al. 2017). In Chapter 2, I showed that summer temperatures were indeed a useful predictor of Arctic grayling distributions at the watershed scale (10s – 100s of km), though data gaps in both thermal habitat use and distributions have been identified at scales useful to planning on-the-ground conservation interventions (< 1 – 10 km; Stamford et al. 2017). Within the Williston Reservoir watershed, these scales represent individual rivers and the discrete reach-scale habitat areas that drive distributions within them. While temperature is often an important factor for riverine fish, other drivers may influence the summer abundance of Arctic grayling in the Parsnip River tributaries. Arctic grayling are territorial and follow multiple instream hierarchies during their summer feeding window. Fish both increase in size with upstream distance (Hughes 1999; Baccante 2010) and form size-based dominance hierarchies around optimal feeding positions in pool habitats (Hughes 1992b). As such, upstream distance and position within the river have been found to be important to predicting the distributions of both Arctic grayling populations (Fitzsimmons and Blackburn 2007; Barker et al. 2011) and other salmonid species (Kaylor et al. 2021; Thompson and Overman 2023). Determining if and how position in the river relates to summer Arctic grayling abundance could help 68 conservation managers prioritize sections of the river where specific work can be conducted to the greatest effect. Defining how position in the river may interact with temperature distributions will provide more fine-scale information on how Arctic grayling use their habitats during their feeding window. Arctic grayling exploitation of their feeding habitats may be further influenced by the presence of potential competitors (mountain whitefish; Prosopium williamsoni) or predators (bull trout; Salvelinus confluentus) that occur in the watershed. The objective of this study was to quantify how Arctic grayling use of thermal and physical habitats compares across the reach and river scales. My hypotheses for this study are as follows; (i) Arctic grayling abundance is driven by temperature at both the river and reach scales; (ii) variation in the abundance index of Arctic grayling can be partially explained by upstream distance at the river scale (10’s of km), but not at the reach scale (< 2 km); and (iii) pool habitats will be important to Arctic grayling at both scales. To test these, I combine snorkel surveys with drone mapping to compare the physical and thermal habitat use of Arctic grayling at the reach scale (< 2 km) and river scale (1 – 52 km). I compare the findings from this chapter with the watershed-scale inferences from Chapter 2 to examine the scale-dependence of temperature as a predictor of Arctic grayling habitat use. Methods Study Area This study occurs at two scales within the Anzac River. The Anzac River is a major tributary of the Parsnip River and has a course of 78 km from its headwaters at 2,495 m to its confluence with the Parsnip River at 730 m. The river drains a mountainous region of the Hart Ranges along the Rocky Mountain trench and has a catchment area of 939 km2. The upper river is characterized by 69 bedrock canyons with a moderate gradient (1 – 2%) which reduces to low-gradient meanders across the unconfined lower river valley as it nears the Parsnip River confluence. Spring freshet in this system is dramatic; snowmelt causes high flows and turbidity that peak soon after ice out and then by late summer gradually reduce to low and clear conditions that persist into the fall. This study occurred over the mid to late summer as flows were tapering off from freshet and Arctic grayling were occupying their summer feeding habitats. Data collection was conducted across two spatial scales: reach and river (sensu Fausch et al. 2002). The reach-scale surveys occurred in two study reaches (named A and B) nested within the river. I selected these reaches based on a combination of site accessibility and the heterogeneity of available thermal habitats within the reach (i.e. tributary inputs and habitats of varying depths and feature classes). These reaches represented the smallest spatial extent in this study and both snorkel and drone sampling were repeated four times throughout the summer of 2022. The river-scale surveys occurred across the first 52 kilometers of the Anzac River and represented the largest spatial extent of the study. The temporal grain of river-scale data was varied; river-scale snorkel surveys were conducted over three years (2018, 2019, 2021) and the river-scale drone survey was conducted once over three days in 2022. Reach A, the lower of the two study reaches, was a 2-km section of the Anzac River defined from river kilometers (Rkm) 37 to 39 (Figure 3.1 and Figure 3.2). It featured two small tributary inputs, the lower of which flows through regions of clearcut land before entering the Anzac and thus brings warmer water into the mainstem (Figure 3.2A). Habitats in reach A included cobble-boulder pools along a bedrock shelf, a long shallow run over bedrock substrates, a prominent deep corner pool (Figure 3.2B), river braids around instream islands, and multiple series of riffle-pool interfaces (Figure 3.2C). A long section of shallow riffles and micro-pools spanned through the lower parts 70 of the reach (Figure 3.2D) before opening to small shallow pools which tailed out at the end of the study area. The elevation at the center of reach A was 797 m. Reach B, the upper of the two reaches, was a 1.1 km section defined from Rkm 46.5 to 47.6 (Figure 3.3). It featured a major tributary input that supported some Arctic grayling in the lower 100 m (Figure 3.3C), a major bedrock chute and series of rapids (Figure 3.3A and Figure 3.3B), series of cobble riffles and pools (Figure 3.3D), a prominent plunge pool below the chute (Figure 3.3A), and a confluence feeding pool below the tributary input (Figure 3.3C). The elevation at the center of reach B is 856 m. Both reaches included a bridge across the mainstem Anzac River. Wetted widths of both reaches remained approximately the same over the sampling season (relative to freshet and in-season discharge spikes observed during other years in this system), though water levels generally dropped as the season progressed (Supplemental Information SI 3.11 and SI 3.12). Snorkel surveys Arctic grayling were counted at both the reach (2022) and river scales (2018, 2019, and 2021) via snorkel surveys. Biweekly reach-scale surveys were conducted four times throughout the summer of 2022 on August 2, August 15, September 2, and September 15 (Table 3.1). Due to logistical challenges, reach B was not surveyed on August 2nd. Reach-scale surveys were constrained to the same spatial extent as a concurrent radio telemetry study (Chapter 4) and then divided into approximately 100-m intervals (subreaches; Figure 3.2 and Figure 3.3) over which all observable Arctic grayling were counted. To determine the potential influences of sympatric species on Arctic grayling habitat use, all observations of a potential predator, bull trout (Salvelinus confluentus), were counted, and the 71 presence/absence of a potential competitor, mountain whitefish (Prosopium williamsoni), were noted. I did not count mountain whitefish because their high abundance limited the capacity to reliably enumerate them. The surveys were conducted by two people swimming each subreach in succession and then taking an average count. Snorkel counts and presence/absence data were integrated with their assigned polygon layers using the R package sf (Pebesma and Bivand 2023; Pebesma 2018; R Core Development Team 2023) and then joined with reach-scale temperature and habitat layers (see below). Subreaches that could not be snorkelled due to very high (e.g. chute section) or low water conditions were excluded from the analysis. River-scale snorkel data were collected as part of a long-term monitoring study on this population and was provided by collaborators John Hagen and Associates. The survey methods were described in their reports (Hagen and Stamford 2023; Hagen and Gantner 2019). Only years with corresponding temperature data available (2018 – 2021, excluding 2020 during which no snorkel surveys were conducted) were used from the snorkel database. Selected river-scale snorkel surveys conducted across 2018, 2019, and 2021 showed consistently high Arctic grayling densities (> 100) in the reaches between Rkms 30 and 47 (Supplemental Information SI 3.1). For the 2018 survey, which covered more of the river than just the index reaches surveyed in other years, high densities were also found in Rkms 23 to 30 (when offset by the variable lengths of the snorkel survey reaches). Rkm 47, just below the chute barrier, consistently had the largest Arctic grayling counts across all years. Counts were lower in surveys below Rkm 16 and above Rkm 48 (above the chute barrier). 72 Drone surveys To quantify the availability of both the physical and thermal habitats in the Anzac River, drone surveys were conducted four times (bi-weekly) at the reach scale and once at the river scale. Surveys were conducted using a DJI Matrice 200 V2 quadcopter drone (DJI, Shenzhen, China) equipped with a Zenmuse XT2 four-band (one band each of red, blue, and green, for optical imaging and one band of forward-looking infrared/FLIR radiometric sensor for thermal imaging; DJI, Shenzhen, China) imaging payload. Drone surveys of each reach were conducted on the same dates as the biweekly snorkel surveys with one exception; the first drone survey in reach A occurred two days earlier (July 31, 2022). The river-scale survey was conducted over three days between August 11 – 14, 2022 (Figure 3.1, Supplemental Information SI 3.2). Surveys were conducted from a series of vantage points accessible through cutblocks along the Anzac River valley. River reaches within 3.3 km of each vantage point were mapped, which depending on the sinuosity of the survey section produced single-flight orthomosaics ranging from 2 to 14 Rkm in length. Surveys were conducted as close as practicable to solar noon at which time the internal heterogeneity of thermal habitats within the wetted channel were the most exaggerated and visible to the thermal sensor (Casas-Mulet et al. 2020; Dugdale 2016, 2015). The river-scale drone survey produced 52 km of continuous imagery over 13 flights of both the physical habitats and local thermal distributions in the Anzac River (Figure 3.1, Supplemental Information SI 3.2). The pool survey identified 482 pools in the first 52 Rkm of the Anzac River which were aggregated by mean size and density per Rkm (Supplemental Information SI 3.18). Mean pool size had a moderately negative correlation with Rkm/upstream distance (‫ = ݎ‬-0.34; 95% CI -0.60 to -0.03) and pool densities weakly increased with upstream distance (‫ = ݎ‬0.31; 95% CI -0.01 to 0.57). 73 While thermal drone studies are now widely used in terrestrial studies (e.g. Larsen et al. 2023; Kays et al. 2019; Chrétien et al 2016), the use of drones for studies focused on temperature in riverine systems is still developing (Dugdale 2016, 2015, Woodget et al. 2017). When the last comprehensive review of drone-based thermal infrared methods was presented by Dugdale et al. (2019), published studies were increasing but still limited at just 12 per year. Following Dugdale’s recommendations, flights were flown in the 0° nadir position to maximize the emissivity of water, though imagery was still subject to variations in temperature from the effects of turbulence, turbidity, or tannin content on emissivity (Dugdale 2016). Other sources of imagery bias could have come from reflections from riparian land cover or atmospheric conditions at the time of sampling. While near-bank reflections cannot be accounted for, temperature logger calibration can be used to account for atmospheric humidity (< 1 °C deviations from FLIR sensor and temperature logger observations). Reach-scale survey B2 had no calibration coefficient available because of a compromised temperature logger. In this case, the calibration coefficient from the reach A survey was used as the two surveys were conducted just one hour apart on the same day under similar atmospheric and camera runtime conditions (Table 3.1). Temperature covariate For the reach scale, I used images from drone surveys to map temperature across the two study reaches. Thermal and physical images produced by the drone surveys were processed with Pix4Dmapper by Pix4D photogrammetry software (version 4.8.3; Educational License; Pix4D S.A., Prilly, Switzerland). Initial processing was completed separately for both the RGB and FLIR rasters using the Advanced Ag RGB and Advanced Thermal Camera processing templates, respectively. Post-processing georeferencing of the RGB to FLIR layers was done in ArcGIS Pro 74 (version 3.1.2; Advanced License; ESRI, Redlands, CA, U.S.A.; Supplemental Information SI 3.3). A full description of the processing methods is given in the Supplemental Information (SI 3.4). FLIR sensors are non-penetrating and thus produce only the surface temperatures of an object (in this case the top 100 μm of the Anzac River water column; Dugdale 2016). To convert the flighttime temperatures of the FLIR raster used in reach-scale analyses to the temperatures during the snorkel surveys (usually <2 hours difference from the time of the drone flights), HOBO MX2201 temperature loggers (Bourne, MA, U.S.A.) were deployed at representative points across each study reach (Figure 3.2 and Figure 3.3). Conversion values for each FLIR layer were calculated as the mean pixel deviations between the temperature logger observations over the duration of the flight and the corresponding surface temperatures recorded by the FLIR sensor at each location and were applied as a linear transformation. As the river-scale drone survey was conducted across multiple days and as such was subject to spatially and temporally varying insolation, cloud cover, temperatures, and local effects, measured temperature values from riverscape FLIR indices were not directly comparable along the longitudinal length of the river, but rather represented local thermal distributions at the time of imaging. For this reason, riverscape FLIR raster values were not used for temperature modeling in the riverscape analyses. Instead, I used mean daily temperature predictions for 2019 and 2021 at a 1-Rkm resolution derived from a spatial stream network model (O’Connor 2023). There was no prediction for 2018, as the specific stream network model is still under development (B. O’Connor, Chu Cho Environmental, unpublished data). Thus, I produced values for 2018 by adjusting the 2021 spatial stream network model predictions by the mean daily variations in temperature 75 observations by temperature loggers installed at static points between the two years (Supplemental Information SI 3.5). I used 2021 as a baseline because it had a similar temperature profile to 2018. Other covariates I used the RGB data from the drone survey to develop further indices of habitat for Arctic grayling. At the river scale, individual pools were surveyed using the drone to create aggregated metrics of pool habitat quality (mean ranked pool size and pool density) at a 1-Rkm resolution (Supplemental Information SI 3.5). Pools were identified as areas where the water became deeper (e.g. changes in water colour or more obscured substrate visibility) and surface flow patterns showed evidence of slowing relative to adjacent flows (less surface disturbance or upwelling patterns). Pools typically formed in the corners of river bends, behind hydraulic features (e.g. boulders or logs), at or below tributary confluences, and at times where the river channel deepened and widened and midchannel flows and substrates changed. Individual pools were scored for size, with their widths and lengths scaled in relation to channel width, so that pools were weighted relative to their local availability within the riverscape. Five size levels were defined (Figure 3.4): - Size 0: Very small, but larger than micropools which cannot be clearly identified at the imagery resolution. Possibly ephemeral features caused by temporary river features such as isolated logs or exposed large boulders creating collecting pools in their wake. Took up < 10% of the overall reach width. - Size 1: Small pools, approximately equal in length and width. Took up 10 – 25% of the overall reach width. - Size 2: Pools which took up 25 – 100% of the overall reach width but were less than two reach widths in length and/or contained inclusions of other habitat types. 76 - Size 3: Large pools which took up the entire width of the reach and had a length of two reach widths or more. - Size 4: Notably large pools or pool complexes which covered the entire reach for three or more reach widths or which showed major widening of the reach when compared to adjacent channels. To determine how upstream distance may influence Arctic grayling abundance at each scale, river kilometers were assigned in 1-Rkm (at the river scale) and 100-m (at the subreach scale) intervals increasing with upstream distance from the confluence of the Parsnip River. At the reach scale only, a habitat layer of three polygon classes (pools, riffles, and runs) was created by manually defining polygons for each habitat class from drone orthomosaic images (Supplemental Information SI 3.6). Modeling Reach-scale snorkel count data were analysed using a negative binomial generalized linear mixed model (GLMM) using the R package glmmTMB (Brooks et al. 2017, R Core Development Team 2023). These modeled the counts of Arctic grayling in each subreach against the suite of environmental covariates. Reach-scale GLMMs took the form: ‫ܴܩ‬௜,௝ ~ܰ‫ܤ‬൫ߤ௜,௝, ∅௜,௝ ൯, (Eq. 1) ଶ ଶ ln൫ߤ௜,௝ ൯ = ߚ଴ + ߛ௜ + ߚଵ ܾ௜,௝ + ߚଶ ‫ݓ‬௜,௝ + ߚଷ ℎ. ‫ݔ‬௜ + ߚସ ‫ݐ‬௜,௝ + ߚହ ‫ݐ‬௜,௝ + ߚ଺ ℎ௜ ‫ݐ‬௜,௝ + ߚ଻ ℎ௜ ‫ݐ‬௜,௝ + ߚ଼ ݇௜ , (Eq. 2) where the count of Arctic grayling ‫ܴܩ‬௜,௝ in subreach ݅ at time ݆ is assumed to follow a negative binomial distribution with dispersion parameter ∅௜,௝ and expected value ߤ௜,௝ ; ߚ଴ is the intercept and ߛ௜ is a random effect of subreach that is normally distributed with mean zero and standard 77 deviation ߪ estimated from the data; ߚଵ is the coefficient for the count of bull trout (ܾ௜,௝ ), ߚଶ is the coefficient for the presence or absence of mountain whitefish (‫ݓ‬௜,௝ ), ߚଷ is the coefficient for the area of habitat ‫( ݔ‬pool, riffle or run) in square meters (ℎ. ‫ݔ‬௜ ); ߚସ and ߚହ are the coefficients for the ଶ polynomial temperature terms ‫ݐ‬௜,௝ and ‫ݐ‬௜,௝ , respectively (polynomial temperature terms were used because I expected a nonlinear relationship between temperature and abundance that would reflect a thermal niche); coefficients ߚ଺ and ߚ଻ are the interactions between the habitat ‫ ݔ‬area and temperature terms and were included as I expected that the effects of temperature would change with habitat areas due to variable rates of heating over the water’s surface; and ߚ଼ is the coefficient for the upstream distance (Rkm) at each 100 m subreach ݇௜ . To avoid the issues of multicollinearity that can arise from using proportion/sum-constrained variables in regression models (Valle et al. 2024), I first conducted a preliminary analysis comparing global models including only one habitat type (ℎ. ‫ )ݔ‬at a time (i.e. run, riffle, and pool; Table 3.2, ‫ܯ‬ௌா் 1). The best model was selected by Akaike’s Information Criterion with correction for small sample sizes (AICc; Burnham and Anderson 2001) and indicated strongly that pool habitats were the best habitat covariate for explaining the variation in the data. Subsequent analyses were conducted using only the area of pool habitats. Reach-scale Rkm and temperature had low correlation (r = -0.10; VIF = 1.01). River-scale snorkel count data were analysed using a negative binomial generalized linear model (GLM) to model the counts of Arctic grayling against habitat covariates using the R package glmmTMB (Brooks et al. 2017, R Core Development Team 2023) with the form: ‫ܴܩ‬௜,௝ ~ܰ‫ܤ‬൫ߤ௜,௝ , ∅௜,௝ ൯, (Eq. 3) ln൫ߤ௝ ൯ = ln൫݈௝ ൯ + ߚఓ,଴ + ߚఓ,ଵ ‫ݑ‬௝ + βஜ,ଶ ‫ݐ‬௝ + βஜ,ଷ ‫ݐ‬௝ଶ + ߚఓ,ସ ‫ݏ‬௝ + ߚఓ,ହ ݀௝ , (Eq. 4) 78 where the count of Arctic grayling ‫ ܴܩ‬over snorkel survey ݅ at location ݆ follows a negative binomial distribution with dispersion parameter ∅௜,௝ and expected value ߤ௝ ; ln൫݈௝ ൯ is the offset term for the log length of a surveyed stream reach; ߚఓ,଴ is the intercept; ߚఓ,ଵ is the coefficient for the river kilometer of upstream distance ‫ݑ‬௝ ;ߚఓ,ଶ and ߚఓ,ଷ are the coefficients for the polynomial temperature terms; ߚఓ,ସ is the coefficient for the mean pool size in each 1-Rkm observation unit ‫ݏ‬௝ ; and ߚఓ,ହ is the coefficient for the pool density in each observation unit ݀௝ . Because temperature and Rkm were strongly correlated in the 2019 data (-0.87), they were not included in the same global model as represented in Eq. 4. Therefore, candidate models were fit using combinations of all terms with either Rkm or the polynomial temperature terms. Candidate models for the reach and river survey data included different combinations of the predictors (yielding a total of 15 and 10 reach and river models, respectively). I used AICc to identify the most parsimonious model. The coefficients of the models that accounted for 95% of total the AICc weight (i.e. 95% confidence set for the best model as per Burnham and Anderson 2002) were averaged to create predictions of Arctic grayling abundance against habitat covariates at each scale. Spatial autocorrelation was checked using plots of residuals by Rkm, and the adequacy of the models chosen to make predictions was checked by assessment of quantile residual plots using the R package DHARMa (Gelman and Hill 2006; Hartig 2022) and comparison to intercept-only models. Results Reach-scale results Arctic grayling snorkel counts in reach A were the highest during survey A1 on August 2, dropping from 50 to 15 by the second survey (Supplemental Information SI 3.13). As the season progressed, 79 numbers increased in reach A with 34 and 40 Arctic grayling counted on September 2 and 15. While there was no snorkel survey in reach B on the first survey date, the reverse trend was generally true; snorkel counts of Arctic grayling were the highest during the second survey on August 15, with numbers decreasing as the season progressed (Supplemental Information SI 3.14). These inverse trends suggest we may have captured the tail end of the upriver summer migration in early August (with fish moving from reach A towards reach B) and the beginning of the fall migration back towards overwintering habitats in the Parsnip River mainstem and lower Anzac River (Blackman 2002; Chapter 2). Model selection indicated that the interactions between pool area and temperature were important covariates for explaining Arctic grayling abundance at the reach scale, as the interaction was included in eight of the top 15 models comprising 0.95 of the cumulative AICc weights (Table 3.2; ‫ܯ‬ௌா் 2). Model averaged predictions of peak counts of Arctic grayling occurred at different temperatures depending on pool size, with peaks occurring at warmer temperatures in small pools and cooler temperatures in large pools (Figure 3.5A). Arctic grayling counts increased strongly with pool area at cool temperatures, showing weaker responses at medium temperatures and a weak reversal in the relationship with area at warm temperatures (Figure 3.5B). Upstream distance was included in the top model and in seven of the top 15 models comprising 0.95 of the cumulative AICc weights (Table 3.2; ‫ܯ‬ௌா் 2). Both reaches A and B showed a slight increase in Arctic grayling index counts with upstream distance (Figure 3.6). While not included in the top model, bull trout counts and mountain whitefish occurrence were present in seven of the 95% confidence set models, indicating uncertainty about their effects on Arctic grayling counts (Table 3.2; ‫ܯ‬ௌா் 2). Model averaged prediction plots showed that Arctic grayling densities weakly decreased as more bull trout were present in each subreach, and weakly 80 increased with the presence of mountain whitefish (Supplemental Information SI 3.7). As the sample size in the reach-scale models was small (n = 59), there were not enough combinations of species, pool area, and temperature in the data to generate reliable predictions of Arctic grayling counts by these effects. Reach-scale models were assessed for fit using quantile residuals and residual autocorrelation plots (SI 3.8, SI 3.9, SI 3.10), and by comparison against an intercept-only model. River-scale results Model-averaged predicted Arctic grayling counts at the river scale were low in the lower Anzac River with a few patches of higher values around Rkms 8 and 17 that increased with upstream distance. Count predictions were greatest for the upper river with the largest predictions occurring around Rkm 45 and decreasing slightly into the highest reaches (Figure 3.8). The abundance of Arctic grayling at the river scale was best described by two count models (AICw > 0.95, sample size n = 39). The top ranked model included river kilometer and mean pool size. The second ranked model also included these predictors as well as pool density (Table 3.3). Modelaveraged predictions revealed that Rkm was the strongest predictor, with Arctic grayling counts increasing markedly with upstream distance (Figure 3.7A). Mean pool size within each Rkm also had a strong positive effect on Arctic grayling counts (Figure 3.7B). However, the relationship between pool density and Arctic grayling counts was weakly positive (Figure 3.7C). Model parameters, quantile residuals, and residual spatial autocorrelation for the river-scale models are available in the Supplemental Information section (SI 3.15, SI 3.16, SI 3.17). 81 Discussion In this study, the area of pool habitats, temperatures (11.0 – 16.0 °C), and upstream distance predicted the abundance of Arctic grayling at the reach scale (100 m – 2 km). This, in part, reflects the finding from Chapter 2 that temperature (11.1 – 17.1 °C) is important for dictating the distribution of Arctic Grayling at the watershed scale (10s – 100s of km). At the intermediate, river scale (1 – 10s of km), metrics of pool habitat (mean ranked size and density per Rkm) were also important to explaining patterns in counts, though upstream distance (Rkm) superseded temperature as a predictor at this scale. Together, these results suggest a degree of scaledependence in explaining Arctic grayling counts during their summer feeding period in the Anzac River. The consequences of multiscale effects in other systems highlight how inferences identified at one scale, while important, may miss other important drivers of habitat use nested at higher or lower hierarchical orders (Fausch et al. 2022; Torgersen et al. 2022). For example, redd density of bull trout spawning in Montana’s Swan River was strongly associated with groundwater upwelling zones in alluvial valley segments (10 – 15 km), but specific sites for redd construction were associated with localized areas of downwelling that supported the oxygenation of eggs (Baxter and Hauer 2000). Spring-run Chinook salmon in the John Day River basin in Oregon showed multiscale selection of thermal habitats. While river-scale (60 km) temperature suggested that fish could avoid critical temperature by congregating in the cool upper reaches, individuals centred their habitat use in warmer reaches, exploiting pool habitats that provided thermal refuge from local temperatures surpassing their upper lethal threshold (Torgersen et al. 1999). These behaviours also varied among basins, where pool habitats were more important in warm streams, and riffle habitats could be proportionally exploited more often in cool streams (Torgersen et al. 1999). 82 Reach-scale results The results from the reach-scale study were aligned with my expectations that temperatures and pool habitats would be important for predicting the abundance of Arctic grayling. Abundance index predictions peaked between 11 – 16 °C, which was consistent with the watershed-scale findings from Chapter 2 (11.1 – 17.1 °C). Arctic grayling abundance showed the strongest response to large pools that supported cool temperatures. Small pools with warmer temperatures were also used but were associated with lower abundances. It is possible that this was an effect of dominance hierarchies, where larger, more dominant fish were able to maintain more favorable positions in large cool pools and smaller fish were forced to occupy suboptimal habitats in warmer pools (this assumes that cooler habitats are preferred, which is confirmed in Chapter 4). While I did not expect upstream distance to be important at this scale, Rkm was included in the best model and in the 95% confidence set. This finding is likely explainable by the morphology of the study reaches. By chance, the lower sections of both study reaches featured large sections of riffle habitat, which did not support large numbers of fish of any species. In addition, the Rkm covariate may be accounting for large differences at the reach scale, with reach A having more fish and being further upstream than reach B. In addition, comparisons between reach-scale global models including either temperature or Rkm (but not both) indicated that temperature was more important to explaining the variation in the data at the reach scale. Future work could be conducted in study reaches with morphologies selected specifically to test this effect further (e.g. study reaches with evenly-spaced pool habitats along their lengths). Water levels in 2022 were low compared to other recent years, which may have emphasized the importance of deeper pool habitats. Rainfall events in 2022 prevented the extreme drought conditions present in 2018 and 2021 (during which these reaches were functionally dewatered) but 83 were low enough to preclude effective snorkel surveys in parts of the shallower reaches. While it is possible that some habitat use occurred in micro-pools along cobbly riffles, they would have been largely unobserved during our snorkel surveys in which these reaches were often surveyed by walking in < 0.5 m-deep water more than they were by swimming, and the dark coloration of Arctic grayling when viewed from above could make spotting them against the substrate difficult. Habitat inferences in the reach-scale study could also have been biased by how the habitat layers were defined. Automated GIS workflows for riverine habitat delineation are limited with respect to the combination of sensors we had on board the drone. As such, habitat areas used in this chapter were defined manually for each study reach and subjective determinations were made based on drone orthomosaics and site knowledge. Model predictions could have been influenced by decisions about where run or riffle habitats ended and pool habitats began. While pool habitats are associated with feeding for Arctic grayling, terrestrial drift production is greatest in riffle habitats (Naman et al. 2017). Future studies could define habitat layers which include interface habitats between pools and riffles, as expected densities of Arctic grayling would likely be highest at these interfaces. River-scale results In this study, the upstream distance (Rkm) covariate outperformed temperature as a predictor of counts in the global models at the river scale. This was unexpected and contrary to my expectations given the emphasis on temperature identified at both the reach (this chapter) and watershed (Chapter 2) scales. Though unexpected, this finding can be explained from an ecological perspective. Upstream distance has been identified as an important factor in other studies of riverine fish growth. While fish growth was not measured in this study, it may have relevance to 84 Arctic grayling distributions during their summer feeding period during which growth is their goal. In the upper Snoqualmie River in Washington, fish were sampled at both the river and section scales, and it was found that while annual fish growth decreased with upstream distance, the rate of summer growth increased (Thompson and Overman 2023). In juvenile Chinook salmon within two Oregon watersheds, upstream distance was related to lower growth rates in the spring, but this pattern reversed in the summer when individuals occupying higher reaches began showing higher rates of growth (Kaylor et al. 2021). These findings related to upstream distance and growth are interesting considering my focus on Arctic grayling during their summer feeding period. Summer distributions of fluvial Arctic grayling are related to dominance hierarchies in which adult and subadult males and females establish and defend feeding territories (Vascotto 1970; Kratt and Smith 1979). These hierarchies are themselves multiscale; individuals distribute themselves within the river following a largerfish-upstream pattern (Hughes and Reynolds 1994; Baccante 2010), and sub-hierarchies emerge around pool habitats in which larger, more dominant individuals maintain optimal feeding positions with respect to net energy income (Hughes and Dill 1990; Hughes 1992a,b; Hughes 1998). While upstream distance has been well-explored in terms of Arctic grayling size patterns, few studies have explicitly looked at patterns in Arctic grayling abundance along this gradient. In the Upper Little Smokey River in Alberta, abundance of both small and large individuals increased with upstream distance up to a point, after which small individuals disappeared and the abundance of large individuals remained constant (Fitzsimmons and Blackburn 2007; Fitzsimmons et al. 2009). In the Peel Watershed in the Yukon Territory, upstream distance was the best predictor of Arctic grayling presence, but it had a negative, though insignificant, effect on abundance (i.e. 85 upstream habitats were more likely to be occupied by dominant fish, but the overall number of fish trended downwards; Barker et al. 2011). While Rkm performed well as a linear predictor in this study, headwater reaches can often be inaccessible to fish due to high gradients and narrow channel widths (Isaak and Young 2023). In my study, the upstream distance metric (Rkm) represented only the first 52 kilometers of the Anzac River as surveyed by drone. At kilometer 47, there is a bedrock chute that has been historically considered a barrier to fish movements (Blackman and Hunter 2001). While it has been more recently understood that Arctic grayling can traverse the chute (which was observed in both the 2018 river-scale snorkel surveys by Hagen et al. [2019]; and in 2022 during exploratory snorkel surveys after August 18 [unpublished data]), it is likely that it still acts as a filter on population abundance (which is supported by the trends shown in Supplemental Information SI 3.1). Indeed, Hagen et al. (2019) found that the abundance of all riverine species decreased above this point and that habitats above the chute quickly became limited in both flows and deep pool habitats. With respect to Rkm as a predictor in the model, my findings may be in part explainable by the defined extent of my survey range. I consider that Rkm would have likely performed better as a polynomial predictor if the headwater portions of the Anzac River had been included. In other words, my study extent likely ended near the upper distribution limit of Arctic grayling in this system, and in doing so the Rkm covariate explained much of the variation in the data, though predictions made with this metric would not likely hold if extrapolated into the headwater reaches. Further considerations of the selection of the Rkm covariate may be explained by the continuum of biophysical processes which occur along longitudinal river gradients (Vannote 1980). The selection of summer feeding sites at the river scale is a function of prey densities. Indeed, multiscale trophic dynamics can complicate the prediction of fish abundance in rivers (e.g. 86 Thompson and Overman 2023; Kelly et al. 2023). Arctic grayling feed on terrestrial drift (McPhail 2007; Stamford et al. 2017), and their success is related to turbidity, the amount of drift produced by different habitat types and conditions, and their fine-scale positions within pool habitats (Hughes 1992a). Invertebrate drift densities along rivers are themselves multiscale processes related to discharge, water quality, light, and riparian condition (Hoover et al. 2007; Kennedy et al. 2014; Riato et al. 2023). While it was considered for this study, sampling invertebrate drift densities at comparable spatial resolutions was deemed logistically prohibitive after consultation with experts (D. Erasmus, UNBC and C. Cena, Environmental Dynamics Inc., Personal Communications). A useful follow-up to this work would be a study focused on the multiscale relationship between the feeding behaviour of Arctic grayling and the density of invertebrate prey. As in other systems with pronounced elevation gradients, temperature is typically negatively correlated with upstream distance in the Anzac River (-0.97 in 2019; though periods of thermal extremes can disrupt this). Considering this, when combined with pool habitats, my models using Rkm may account for two of the main drivers of Arctic grayling habitat use during summer feeding. Warming along the longitudinal length of rivers is a dynamic process dependent on tributary inputs, air temperatures, velocity, precipitation, and elevation (Fullerton et al. 2015). Water temperatures can further vary at different spatial resolutions, and in hot years during which asymptotic warming can disrupt otherwise linear temperature profiles, fine-scale use of thermal refuge habitats can drive the distribution of species (Torgersen et al. 1999; Fullerton et al. 2015; Fullerton et al. 2018). My results have important implications for conservation planning. If upstream distance and pool habitats are indeed the best explanatory variables for explaining summertime river-scale Arctic grayling abundance in the Anzac River, then this implies that Arctic grayling abundance may be 87 driven more by their position within the river than water temperature (bearing in mind that upstream distance does have some effects of temperature gradient). Indeed, Arctic grayling distributions at the river scale are driven by dominance hierarchies for optimal feeding habitats (Hughes and Dill 1990). While some inter-annual redistribution of individuals within rivers is necessary to maintain population size gradients as individuals grow (Hughes 1992a; Christie et al. 2010), site fidelity is generally high in fluvial Arctic grayling populations. In a radio tagging study, Arctic grayling were often observed moving back to the same sites at which they were tagged during the previous summer (Blackman 2002), a pattern which was reinforced in exploratory analyses of the acoustic telemetry dataset used in Chapter 2 (Martins et al. 2022). If site fidelity and maintaining these hierarchies at the river scale are more important drivers of their summer abundance than temperature, then conservation planning based on temperature control alone may miss important aspects of Arctic grayling conservation in this system. I emphasize that my results do not show that temperature is unimportant at the river scale. My results only demonstrate that Rkm was a stronger predictor of river-scale abundance in the Anzac River in 2018, 2019, and 2021. In this study, I used modelled temperatures in the river-scale study (O’Connor 2023), but the temperature models for 2018 (during which the broadest extents of the river-scale snorkel surveys were conducted) are still under development (B. O’Connor, Chu Cho Environmental). The summer of 2018 was warm, and data from instream temperature loggers most closely matched those observed in 2021; as such the 2018 temperature data used in this model are an adjusted representation of the 2021 model. In the 2021 model, the correlation between temperatures and Rkm are lower than the one found in 2019 (a cooler and more ‘typical’ year; Supplemental Information SI 3.19). For this reason, temperature outliers in the 2021 data, and their 88 reflections in the 2018 data, may have masked the importance of temperature in the river-scale models. Conservation planning in rivers often focuses on the coolest summer habitats, and there is a need to evaluate the contribution of warmer parts of the river to population success (Kaylor et al. 2021). As in the Chinook salmon study by Torgersen et al. (1999), the scale-dependence of temperature is likely important in this system. Reach-scale findings suggest that local exploitation of temperatures is related to the area of pool habitats within the reach. As the importance of pools was maintained at both the river and reach scales, it may be that river-scale habitat use as related to the local availability of structured habitats, which support both thermal refuge and foraging opportunities. As such, conservation interventions in this system (e.g. land procurement for restoration or protection) should, in addition to the important conservation of headwater reaches which support cold water inputs, also consider reaches within the river which support highly structured instream habitats. Given my findings, the combined prioritization of both habitats would lead to better outcomes for Arctic grayling conservation in the region. In this chapter, I identified a degree of scale dependence in predicting the abundance of Arctic grayling during their summer feeding window. As supported by other studies on this species (e.g. Zemlak and Langston 1998; Blackman et al. 2002), pool habitats were found to be important at both the reach and river scales. Temperature, a potential limiting factor in many freshwater species (Dudgeon et al. 2006), was found to be important at the reach scale presented in this chapter and the watershed scale presented in Chapter 2. However, at the river scale, upstream distance superseded temperature as an environmental predictor, a novel finding that highlights the importance of multiscale investigations in conservation planning. 89 In addition to the conservation of headwaters, this work highlights another tool conservation planners can use to improve the outcomes for this species. In addition to the protection of coldwater habitats, Rkm may be useful when considering the conservation of Arctic grayling in the Williston Reservoir at large by highlighting candidate areas for spatial fisheries closures during the summer season. The use of such a simple metric, that does not require in-depth characterization of the river habitat, would be advantageous given that the system is largely inaccessible. Indeed, Stamford et al. (2017) called for a GIS-based approach to Arctic grayling habitat studies that could benefit from this finding. However, given that the Parsnip River watershed and Anzac River are merely one catchment of a system covering thousands of square kilometers across northern British Columbia, the suitability of this metric would have to be explored further. In an era of increasingly uncertain climate stability, the temperature related findings of this study broadly reinforces that the conservation of headwaters is important in systems which support Arctic grayling. However, site fidelity and territoriality may drive the summer distributions of this population and could increase their thermal risks as they defend territories in shrinking areas of thermal refuge. Consequently, the potential effects of the distribution of thermal refuges within the riverscape should be studied and are likely also important to making accurate predictions within extreme thermal years. 90 Reach A B A B A B A B Survey 1 1 2 2 3 3 4 4 2022-09-15 12:28 - 12:38 2022-09-15 11:33 - 11:41 2022-08-15 15:14 - 15:23 2022-08-15 14:27 - 14:34 2022-09-02 12:36 - 12:44 2022-09-02 11:46 - 11:56 2022-07-31 14:05 - 14:16 2022-07-31 13:04 - 13:19 Flight time 2022-09-15 13:05 - 14:00 2022-09-15 14:45 - 15:33 2022-08-15 12:00 - 13:05 2022-08-15 16:10 - 17:00 2022-09-02 13:30 - 14:20 2022-09-02 15:00 - 15:55 2022-08-02 11:00 - 12:30 No snorkel survey Snorkel time 9.7 9.9 14.7 15.9 11.4 11.2 14.9 14.9 FLIR temp (°C) 9.7 9.8 13 16 12.1 14.4 12.3 - HOBO temp (°C) 91 0 -0.1 -1.7 0.1 0.7 3.2 -2.6 - Calibration (°C) Table 3.1. Table of each drone and snorkel survey conducted in the reach-scale study within the Anzac River, as well as the calibration values between the FLIR sensor and HOBO temperature loggers at the time of the flight. Chapter 3 Tables . Table 3.2 AICc table for the reach-scale modeling study, where GR is the count of Arctic grayling in each 100-m subreach, ‫ ݐ‬+ ‫ ݐ‬ଶ is the polynomial temperature term, ℎ. ‫ ݔ‬is the area in ݉ଶ of each habitat type ℎ. ‫( ݔ‬pool, riffle or run) within each subreach, ܾ is the count of bull trout during each survey, ‫ ݓ‬is a binary presence/absence parameter of mountain whitefish occurrence in each subreach during each survey. All models include a random effect of subreach on the intercept. Presented are the number of parameters in the model (K), the log likelihood (logLik), the AIC score with small sample size correction (AICc), the difference in AICc from the top model (dAICc), and the model weights (AICcwt). Models which are bolded in ‫ܯ‬ௌா் 2 were included in the 95% confidence set. MSET 1: Select global models for best habitat type Model k logLik AICc dAICc AICcwt GR ~ b + w + h.pool x (t + t2) 10 -144.96 314.49 0.00 1.00 2 GR ~ b + w + h.riffle x (t + t ) 10 -151.11 326.81 12.31 0.00 GR ~ b + w + h.run x (t + t2) 10 -151.31 327.21 12.72 0.00 MSET 2: Select candidate models for pool habitats Model k logLik AICc dAICc AICcwt GR ~ h.pool x (t + t 2) + k 9 -144.12 309.90 0.00 0.24 GR ~ h.pool x (t + t 2) 8 -145.85 310.59 0.68 0.17 2 GR ~ w + h.pool x (t + t ) + k 10 -143.26 311.10 1.19 0.13 GR ~ w + h.pool x (t + t 2) 9 -145.14 311.95 2.04 0.09 GR ~ b + h.pool x (t + t ) + k 10 -144.11 312.79 2.89 0.06 GR ~ h.pool + k 5 -150.88 312.89 2.98 0.05 GR ~ k 4 -152.15 313.04 3.14 0.05 6 -149.76 313.13 3.23 0.05 9 -145.85 313.38 3.47 0.04 11 -143.17 313.97 4.06 0.03 GR ~ b + w + h.pool x (t + t ) 10 -144.96 314.49 4.59 0.02 GR ~ b + h.pool + k 6 -150.76 315.14 5.24 0.02 GR ~ h.pool 4 -153.40 315.54 5.64 0.01 GR ~ b + w + h.pool + k 7 -149.76 315.71 5.81 0.01 GR ~ 1 2 -156.43 317.07 7.17 0.01 GR ~ (t + t2) 5 -153.08 317.29 7.38 0.01 2 GR ~ w + h.pool + k 2 GR ~ b + h.pool x (t + t ) GR ~ b + w + h.pool x (t + t 2) + k 2 92 Table 3.3. AICc table for the river-scale modeling study, where ݈ is the offset logged length of each snorkel survey in km to account for the variable length (1.2 – 5.4 Rkm) of the snorkel surveys (݈ is included in all models), ‫ ݑ‬is the upstream distance from the confluence of the Anzac River with the Parsnip River in km, ‫ ݐ‬+ ‫ ݐ‬ଶ is the polynomial temperature term, ‫ ݏ‬is the mean size of all pools within each Rkm, and ݀ is the density of pools within each Rkm. Presented are the number of parameters in the model (K), the log likelihood (logLik), the AIC score with small sample size correction (AICc), the difference in AICc from the top model (dAICc), and the model weights (AICcwt). Models in bold were included in the 95% confidence set. Model k logLik AICc dAICc AICcwt GR ~ l + s + u 4 -198.29 405.76 0.00 0.77 GR ~ l + s + u + d 5 3 4 -198.20 -204.12 -204.11 408.22 414.92 417.39 2.46 9.16 11.63 0.22 0.01 0.00 4 2 5 5 -210.58 -213.84 -210.54 -210.57 430.33 432.01 432.89 432.95 24.57 26.25 27.13 27.19 0.00 0.00 0.00 0.00 GR ~ l + d 6 3 -210.48 -214.71 435.58 436.12 29.82 30.36 0.00 0.00 GR ~ l + s 3 -217.36 441.41 35.65 0.00 GR ~ l + u GR ~ l + u + d GR ~ l + (t + t2) GR ~ 1 GR ~ l + (t + t2) + s GR ~ l + (t + t2) + d GR ~ l + (t + t2) + d + s 93 Chapter 3 Figures 94 Figure 3.1. Orthomosaic image of the Anzac River from river kilometers (Rkm) 0 - 51. This figure is broken into three segments for the lower (left), middle (center), and upper Anzac River (right). River segments have been rotated per the corresponding compass orientations to fit in one image. Purple regions highlight the two nested study reaches in this study. 95 Figure 3.2. Orthomosaic image of reach “A”, the lower of the two reaches from the reach-scale study from Rkm 37 39. It features two small tributary inputs, the lower of which enters the Anzac at temperatures warmer than mainstem temperatures after flowing through regions of clearcut land use (Callout A). Habitats in reach A were varied, including cobble-boulder pools along a bedrock shelf, a long shallow run over bedrock substrates, a prominent deep corner pool (Callout B), river braids around instream islands, and multiple series of riffle-pool interfaces (Callout C). A long section of shallow riffles and micro-pools spans through the lower parts of the reach (Callout D) before opening to small shallow pools and tailing out of the study area. Subreaches from the snorkel surveys (e.g. AR01) are labeled by the black outlines with white text. Imagery has been slightly rotated per the compass rose from true north to fit on the page. 96 97 Figure 3.3. Orthomosaic image of reach “B”, the upper of the two reaches, was a 1.1 km reach defined from Rkm 46.5 - 47.6. It featured a major tributary input which supported some Arctic grayling use in the lower 100 m (Callout C), a major bedrock chute and series of rapids (Callouts A, B), series of cobble riffles and pools (Callout D), a prominent plunge pool below the chute (Panel A), and a prominent feeding pool below the tributary confluence (Callout C). Subreaches from the snorkel surveys (e.g. BR01) are labeled by the black outlines with white text. 98 99 100 Figure 3.4. The five size classes of pool habitats defined at the river scale. Size 0: Very small, but larger than micropools which cannot be clearly identified at the imagery resolution. Possibly ephemeral features caused by temporary river features such as isolated logs or exposed large boulders creating collecting pools in their wake. Took up < 10% of the overall reach width; Size 1: Small pools, approximately equal in length and width, taking up 10 - 25% of the overall reach width; Size 2: Pools which took up 25 - 100% of the overall reach width but were less than two reach widths in length and/or contained inclusions of other habitat types; Size 3: Large pools which took up the entire width of the reach and had a length of two reach widths or more; Size 4: Notably large pools or pool complexes which covered the entire reach for three or more reach widths or which showed major widening of the reach when compared to adjacent channels. 101 Figure 3.5. Model averaged predictions of the effects of temperature (Panel A) and pool size (Panel B) on Arctic grayling abundance at the reach scale. An interaction effect between these two predictors was found, with a peak response in abundance at temperatures of 12.5 °C in large pool habitats, 13.5 °C in medium pool habitats, and 15.3 °C in small pool habitats (Panel B). The response in abundance to pool areas was highest at cool temperatures, weakly increased with pool area at medium temperatures, and weakly declined with pool area at high temperatures. The three panels of each series of plots represent the 0.25, 0.50, and 0.75 quantiles of temperatures and pool sizes, respectively. Plot points indicate the observed data in each modeled combination. Shaded regions represent the 95% CI for the predictions. 102 Figure 3.6. Model averaged predictions of Arctic grayling abundance as a function of river kilometer from the reachscale study. Reach B showed a stronger response with higher Arctic grayling counts than Reach A. Plot points indicate the observed data from the snorkel surveys. 103 Figure 3.7. Model-averaged predictions of Arctic grayling abundance as a function of river kilometer (Rkm; Panel A), mean pool size per Rkm (Panel B), and pool density per Rkm (Panel C) from the river-scale study. The strongest predictor was Rkm with mean pool size showing a similarly strong effect, whereas the effect of pool density was small. 104 105 Figure 3.8. Model averaged predictions of Arctic grayling abundance in the Anzac River at a 1 Rkm resolution. The two models used for predictions included terms for upstream distance (Rkm), mean pool size per Rkm (top and second ranked models), and pool density per Rkm (second ranked model only). Trends in predicted Arctic grayling abundance were low in the lower Anzac River with a few patches of higher densities around Rkms 8 and 17 and increased with upstream distance. Abundance predictions were highest in the upper river with the highest predictions occurring around Rkm 45 and decreasing slightly into the highest reaches. Chapter 4 - Arctic grayling (Thymallus arcticus) show diel and seasonal patterns of behavioural thermoregulation Abstract 1. As ectotherms, freshwater fishes must balance their thermal preferences against other life history goals (e.g. foraging) that require spending time in waters beyond their preferred thermal range. Achieving this balance requires fish to exploit spatiotemporal heterogeneity in their thermal habitats. Behavioural thermoregulation is critical for freshwater fishes to buffer the immediate effects of a warming climate, though few studies have quantified the phenomenon in free ranging fish. I used a population of Arctic grayling (Thymallus arcticus) to test how the spatiotemporal distribution of heterogeneous thermal habitats related to effective behavioural thermoregulation. 2. I used shuttlebox experiments to determine baseline information about both the thermal preference and rate of heat transfer in Arctic grayling that I compared against the behaviours of tagged, free-ranging individuals. Through the first experiment I determined that the thermal preference range for Arctic grayling in the Anzac River ranged from 10.113.0 °C (12.0-14.5 °C for males and 8.8-12.0 °C for females). The second experiment found that the rate of heat transfer was influenced by the condition factor of the fish and whether the body temperature was lower or higher than ambient water temperatures. 3. I quantified behavioural thermoregulation of the free-ranging individuals by calculating metrics of how body temperatures compared to both the thermal preference of the individual and the availability of the thermal preference range in the environment. Using regression analyses, I found that Arctic grayling used behavioural thermoregulation to 106 maintain their body temperatures around their preferred range, and these behaviours were more effective in habitats that were warmer than their thermal preference range. Exploitation of preferred thermal habitats quickly decreased as waters warmed, providing evidence that diel trade-offs with thermal preference occured in this population. 4. Generalized additive mixed models indicated that the effectiveness of behavioural thermoregulation was variable over the 24-hour cycle. The most effective use of this strategy occurred during summer afternoons when temperatures were highest. This pattern reversed at the beginning of the fall when temperatures started to drop (~ September 4). Small energy expenditures increased thermoregulatory effectiveness, but there was a diminishing return where continued expenditures did not result in more effective thermoregulation. Model selection suggested that body condition and the patchiness of thermal habitats may be important to this strategy, though these effects were uncertain. Models with condition factor as a covariate suggested that fish with lower condition factors, which have less thermal inertia than fish with higher condition factors, used behavioural thermoregulation more effectively. Behavioural thermoregulation was marginally more effective in patchier thermal habitats. Results showed that Arctic grayling may be among a subset of species which thermoregulate in only one direction (cooling). Introduction Temperature mediates the ecophysiology of freshwater fishes: it can act as a cue for migration (Elsner & Shrimpton 2019), spawning (Hubert et al. 1985; Shuter et al. 2012), and diel foraging patterns (Armstrong et al. 2013), dictate the times conducive to (and intensity of) activity (Gunderson & Leal 2016; Abram et al. 2017), and influence survival and reproductive success (Martins et al. 2011; Dahlke et al. 2020). As freshwater fishes are ectothermic, their body 107 temperatures (and all subsequent physiological functions) are directly determined by their position within the thermal environment. At the population and species levels, adaptations to changing thermal environments occur over years to decades through selective pressures on their underlying genetics (Kearney et al. 2009; Kelly 2019). In contrast, individuals cope with changes in their thermal environment in the intermediate-term (days to months) through physiological acclimation (Sears et al. 2019). Over the short-term (minutes to hours), individuals must rely on behavioral thermoregulation (i.e. seeking out cooler or warmer habitat) to regulate their body temperatures around a thermal preference range and avoid unfavourable critical temperatures (Angilletta 2009; Farrell et al. 2008). Behavioural thermoregulation in most ectotherms can only be effective in spatially heterogeneous thermal habitats in which a range of temperatures are available for individuals to exploit (Sears et al. 2016; Sears and Angiletta 2019). The effectiveness of behavioural thermoregulation is quantified by measuring how both body and ambient temperatures deviate from an established thermal preference range (Hertz et al. 1993; Blouin-Demers and Weatherhead 2001; BlouinDemers and Nadeau 2005). To-date, most investigations into the effectiveness of behavioural thermoregulation have been assessed on terrestrial ectotherms in which both the spatial structure and heterogeneity of temperatures at localized spatial scales were important to the effectiveness of these behaviours (e.g. Vickers et al. 2011; Sears and Angiletta 2015). In freshwater habitats, temperatures are not as varied as in terrestrial habitats, but they can also vary across localized spatial scales more strongly than previously thought (Kurylyk et al. 2015; Dzara et al. 2019). For example, it is well known that temperature along the water column of temperate lakes during the summer can range widely from about 20°C at the surface to about 4°C in deep waters (Boehrer and Schultze 2008). Less appreciated is the spatial variation in water 108 temperatures that occur along rivers and across their channels. Riverine fishes can encounter water temperature differences ranging from a few (1 – 2 °C) to several (3 – 6 °C) degrees over 10 – 100s of metres due to the input of coldwater tributaries and groundwater exchange (Pincebourde et al. 2016; Kurylyk et al. 2015; Fausch et al. 2002). Water temperature in freshwater habitats can also change markedly over short (hourly, daily) and long (within and between seasons) temporal scales (Caissie 2006), providing fish with a spatiotemporally complex distribution of heterogeneous thermal habitats in which to regulate their body temperatures. As daily and seasonal temperature cycles progress, the associated changes in the availability of thermal habitats can influence how thermoregulatory behaviours are exploited over time (Karameta et al. 2023). Theory suggests that the energetic costs of thermoregulation increase in structurally clumped (vs. dispersed) thermal habitats with high variability in temperatures (Sears et al. 2016). However, allocation of energy resources towards thermoregulation may also depend on how suitable thermal habitats overlap with areas suited to other goals (e.g., foraging, reproduction, or predator avoidance; Veech 2021). For example, juvenile broad-headed snakes (Hoplocephalus bungaroides) showed a 57% reduction in time spent within preferred body temperature ranges when exploiting habitats that reduced the risk of predation (Webb and Whiting 2005). While studies have long acknowledged the fitness costs associated with the tradeoffs in behavioural thermoregulation (e.g. Huey and Slatkin 1976), assessing the energetic costs of this strategy has received relatively less attention. This may be attributed in part to mixed results between field and laboratory studies on thermoregulation (Basson et al. 2017). The optimality model of thermoregulation (Huey and Slatkin 1976) suggests that thermoregulatory behaviours will be energetically costly to maintain and that individuals in favorable thermal habitats will have to expend less energy. However, studies in terrestrial ecosystems have not always found this to be the 109 case. Experimental studies found behavioural thermoregulation incurred low energetic costs in turtle embryos (Zhao et al. 2014). Studies on lizards found that the same body temperatures were maintained regardless of habitat quality and individuals invested heavily in thermoregulation even in scenarios where thermoconformity would have resulted in significant energy savings (Basson et al. 2017). As individuals navigate the tradeoffs presented by their environments, the duration of excursions beyond suitable thermal habitats into habitats suited for other goals can vary with body size as a result of thermal inertia, where the body temperature of larger individuals takes longer to change and equilibrate to a different ambient temperature (Iverson et al. 1997; Stevenson 1985; Pépino et al. 2015). Thermal inertia is closely linked to the daily thermal experience of small ectotherms and can drive the concept of ‘inertial homeothermy’ in very large ectotherms (Stevenson 1985). In turtles, the effect is behind both diel and seasonal variation in nesting timing, where larger individuals can retain their body heat from the aquatic environment for longer when making excursions into cooler air temperatures to lay eggs (Connoy et al. 2020). Thermal inertia is the physical property that can contribute to emergent patterns in thermal ecology related to body size and temperature. In northern water snakes (Nerodia sipedon), body size and temperature play a predominant role in reproductive success within micro and macro scales (Brown and Weatherhead 2000). There is evidence, therefore, that body size and temperatures are linked to distribution patterns in reptiles (Angiletta et al. 2004), though this effect has not been well explored in fishes (Amat‐Trigo et al. 2023). The heat transfer rate of an organism, or how quickly internal body temperatures change with exposure to variable ambient temperatures (Huey and Slatkin 1976; Blouin-Demers and Nadeau 2015; Wolkovich et al. 2014; Pépino et al. 2015), is also linked to body size and thermal inertia. 110 While body size can affect the rate of heat transfer in fishes (Steven and Fry 1974; Fechhelm and Neill 1982; Weller et al. 1984; Pépino et al. 2015), they also transfer heat through direct convection at the gills (Stevens and Sutterlin 1976). The proportion of heat transfer through gill and body can vary by species (Stevens and Sutterlin 1976; Masser and Neill 1986). Larger fish with lower heat transfer and higher thermal inertia can take longer excursions into unsuitable thermal habitats for other activities (e.g. foraging), which is supported by simulation experiments (Pépino et al. 2015). In this chapter, I quantifed how effectively adult Arctic grayling (Thymallus arcticus) use behavioural thermoregulation to maintain their body temperatures around their thermal preference range. To do so, I studied a population in the Anzac River, within the Parsnip River watershed of north-central B.C., Canada, during their summer feeding period, when water temperatures are at their highest (Stamford et al. 2017). I tested the hypotheses that (1) thermoregulation was more effective when the riverine habitat was more thermally heterogeneous; (2) effective behavioural thermoregulation required an increased level of activity (i.e. more energy expenditures); and (3) the effectiveness of thermoregulation varied with the pronounced daily temperature cycles in rivers. To do this, I first determined a summer thermal preference range for adult Arctic grayling caught in two reaches of the Anzac River. To quantify how effectively they use behavioural thermoregulation during hot summer months, I then related this thermal preference range to the body temperatures of free-ranging Arctic grayling in the Anzac River. Methods Study Area The methods in this chapter are divided into two categories: methods for the ex-situ laboratory studies and methods for the in-situ study on free-ranging Arctic grayling. Laboratory studies were 111 conducted in a field laboratory set up by Goose Lake near the confluence of the Anzac and Parsnip Rivers. The field laboratory was set up in a large wall tent powered continuously by two 4500w generators hooked in a parallel circuit. In-situ studies were conducted in two study reaches within the Anzac River (named reach A and B) as described in Chapter 2. Fish capture Arctic grayling were captured by angling in the Anzac River using barbless dry flies and sampled for sex, fork length (mm), weight (g), and dorsal fin length (mm). Sex was determined visually based on morphometrics (prominent dorsal fins and broad head indicating males) when possible. For individuals where sex could not be clearly determined, sex was determined by the ratio of dorsal fin length to fork length – individuals were estimated as males if the dorsal length exceeded 50% of fork length. While the ratio of dorsal length and body length has not been fully evaluated as a method for sex determination, a recent study showed that posterior dorsal height can be used for sex determination (Samuel et al. 2023). After capture, condition factor of each fish was calculated based on their length and weight as ‫ܨ‬௜ = ଵ଴రௐ೔ ௅య೔ , (Equation 4.1) where ܹ௜ and ‫ܮ‬௜ are respectively the weight (g) and the fork length (mm) of fish ݅ (Williams 2000). Fish used in the shuttlebox experiments were transported in an aerated cooler to the field laboratory. Capture, handling and tagging (for both the shuttlebox experiments and the radio tagging study) were authorized under the BC Ministry of Forests scientific fish collection permit no. PG22-738069 and done in accordance with protocol no. 2021-05 approved by the University of Northern British Columbia’s Animal Care and Use Committee. 112 Shuttlebox experiments Shuttlebox studies are laboratory-based experimental approaches to establish a range of environmental preferences for an individual while removing the potential effects of factors such as predators, prey, and forage availability (Christensen et al. 2021). The shuttlebox system used in this study (Loligo Systems, Viborg, Denmark) consisted of two round tanks (each 1 m diameter) connected by a passage chamber (20 cm length ´ 15 cm width) to allow individuals to move freely between the tanks. For thermal preference experiments, the shuttlebox system maintained a temperature differential (e.g. 2°C) such that fish could exploit a warm or a cool tank. The system’s software (Shuttlesoft, Loligo Systems, Viborg, Denmark), which was connected to cameras (that tracked the location of the fish between tanks), thermometers, water pumps, hoses, coils, and external heating and cooling baths, automatically warmed the water in both tanks if the fish was in the warm tank and cooled the water in both tanks if the fish was in the cool tank. The fish was therefore able to control the temperature of its body by moving between tanks, revealing its thermal preference over the duration of the experiment (typically 24h or more for each fish). I conducted this thermal preference experiment on adult Arctic grayling caught during their summer feeding period in the Anzac River. Before the experiments, the shuttlebox was disinfected by flushing the pumps and tanks with a solution of Virkon Aquatic disinfectant (Lanxess, Cologne, Germany) for 15 minutes. A 15% water change was done between successive experiments. All fish were released back at their capture locations after the experiments concluded. 113 Thermal preference experiments For each of the thermal preference experiments, I first randomly assigned a tank as either warm or cool. I then set their temperatures to ± 1 °C, respectively, from ambient river temperature at the time and location a fish was captured to create a temperature differential of 2°C between tanks. After setting the system to change temperature at a maximum rate of 2°C / hour, I placed a fish randomly in either the warm or cool tank for a 24-h thermal preference trial. Data from the first four hours of each thermal preference experiment were removed as an acclimation period while the individual explored and learned the shuttlebox environment (Christensen et al. 2021). The remaining 20 hours of continuous temperature data were used to compute an individual’s thermal preference (TPREF) as the 50th quantile (i.e. median) of occupied temperatures; and its preferred temperature range (TSET) as the 25th and 75th quantiles of the occupied temperatures (Christensen et al. 2021). A total of seventeen individuals were used in thermal preference experiments, but data for seven of these fish were not used in this study as they failed to show evidence of learning the shuttlebox system (i.e. settling in one tank for 20 hours, artificially driving thermal preference estimation in one direction). Of the ten successfully trialed, six were females and four were males. Heat transfer experiments For the heat transfer experiments, I blocked the central passage and manually adjusted the temperature of each shuttlebox tank so that each fish would be exposed to a pre-defined shift in temperature. I then monitored the fish for 45 minutes while their body temperatures converged to the water temperature (Table 4.1). At the start of each experiment, individuals were randomly assigned to either the warm tank or the cool tank. The starting tank was set to 2 °C away from ambient river temperatures in the respective direction. Every 45 minutes, the fish would be shuttled to the opposite tank (i.e. cooling experiments were immediately followed by warming experiments 114 at the differential intervals listed below; Table 4.1). Temperature differentials were progressively narrowed as the experiments were conducted (differential intervals for the five trials in each experiment were 2.0, 2.0, 1.5, 1.0, and 1.0 °C; Table 4.1). While a heat transfer experiment was being conducted in one tank, the shuttlebox system would adjust the temperature in the opposing tank to the starting temperature for the subsequent experiment. Before being placed in the shuttlebox tank for the heat transfer experiments, Arctic grayling were equipped both internally and externally with clock-synchronized temperature loggers (Star ODDI DST Nano-T, Garðabær, Iceland) which recorded body and ambient water temperatures every five seconds over each 3.75 h experiment. Internal loggers were attached to a spaghetti tag, inserted esophageally using a plastic applicator coated in Vaseline, and removed after the experiment (Supplemental Information SI 4.5). External loggers were mounted on an anchor tag (FLOY Mfg. T-Bar Anchor, Seattle, U.S.A.), which was inserted under the dorsal fin. The heat transfer experiments were conducted immediately following the thermal preference experiments for seven of the ten fish. Three additional experiments were done using purpose-caught fish after three of the thermal preference fish expelled their internal or external tags midway through the experiment. Heat transfer experiments were conducted after the thermal preference experiments rather than concurrently with them as it was observed during early experiments (data not used) that the loggers affected behaviour in the thermal preference experiments (e.g. fish would stop exploring the tanks or try to remove the esophageal logger using the corner of the shuttlebox tanks and the passage, confounding any results regarding thermal preference, which assume that behaviour is based on selection among available thermal environments alone). 115 Biologging with radiotelemetry Adult Arctic grayling were tagged river side with radio transmitters (model MCFT3-SO-TPA-L, Lotek Wireless). The tags were programmed to alternately transmit a measurement of temperature, pressure (not used in this study) and activity every 5 seconds, producing a 15-second transmission cycle for each metric. The activity sensor measured the rate of change in body acceleration (i.e. jerk; Lennox et al. 2019) over the 15-second sampling interval. The rate of change in body acceleration is a measure of activity that can be used as a proxy for the energetic costs of thermoregulatory behaviours, with higher values of activity assumed to indicate more energy expenditure (Lennox et al. 2023). Radio tags were affixed to the body at the base of the dorsal fin using an interrupted suture technique using monofilament suture line against thin plastic plates spaced from the body with rubber foam padding (Crook 2004). To monitor the deployed radio tags, an array of solar-powered radio receivers (model SRX-1200D2, Lotek Wireless) was installed along each study reach at vantage points on top of bluffs and cliffs to provide maximum line-of-sight coverage for the two (upstream and downstream-facing) receiver antennae (Figure 4.1). Radio tag data recorded by the radio receivers were indexed by study hour with hour zero starting at 2022-07-31 00:00:00 PDT and hour 1,120 ending after 202209-15 16:00 PDT. Between the two study reaches, 50 adult Arctic grayling were tagged with radio transmitters (Supplemental Information SI 4.6). Data were filtered to remove low-power detections (RSSI < 125) that were mostly due to noise caused by the signal echoes and collisions in the bedrock morphology of the study reaches. Further noise may have been present from high radio activity from industrial activity in the area. Tags which produced very limited data after filtering (< 10 detections over the study period with no sequential detection series occurring within 5 116 minutes of each other) were removed, and remaining tags were filtered with a fine-scale filter which removed isolated detections outside of the core detection series. Quantification of thermal habitat The distributions of water temperatures in the study reaches were measured using a drone equipped with a forward-looking infrared (FLIR) camera (as described in Chapter 2). To both validate the thermal rasters collected by the drone and measure temporal temperature trends over the study period, a temperature logger array was deployed in the study reaches (Figure 4.1). As the effectiveness of thermoregulation depends on the availability and the distribution of thermal habitats within the river (Shepard et al. 2013; Sears et al. 2015), the patchiness of thermal habitats was defined using hourly interpolated FLIR rasters which were created using an hourly adjustment value of mean temperatures from temperature loggers in the river with the R package terra (Hijmans 2023; Supplemental Information SI 4.2). The hourly metric of the patchiness of temperatures from these rasters was computed using the R package landscapemetrics (Hesselbarth et al. 2019) as ܲ= ே ஺ × 10଺ (Equation 4.2) where ܲ is the metric of patchiness, ܰ is the number of patches (as determined by contiguous temperature patches in which a temperature pixel is connected to similar temperatures in eight directions), and ‫ ܣ‬is the wetted area of each reach (Hesselbarth et al. 2019). ܲ increases as the thermal landscape becomes patchier. 117 Data analyses Heat transfer coefficient Data from the heat transfer experiments were used to estimate the heat transfer coefficient, k, for adult Arctic grayling. Coefficient k is a metric measuring the instantaneous rate of change in body temperature based on the difference between body temperature (ܶ஻ ) and the ambient water temperature (ܶ஺ ) (Pépino et al. 2015). The coefficient is a parameter in the following heat transfer equation based on Newton’s Law of Cooling: ் ܶ஻೔,೟ = ൬ܶ஻೔,೟షభ − ܶ஺ ௜,௧ିଵ − ௞೘ ൰ ೔,೟ ି௞೔,೟ ୼௧ ்௠ + ܶ஺ ௜,௧ିଵ + ௞ + ߝ௜,௧ , ೔,೟ (Equation 4.3) where ܶ஻೔,೟ and ܶ஻೔,೟షభ are the body temperature of individual ݅ (from the internal tag) at times ‫ݐ‬ and ‫ ݐ‬− 1, respectively; Δt is each 5-second data logging interval, ݇௜,௧ (units of min–1) is the coefficient of heat transfer of individual i for the interval between ‫ ݐ‬− 1 and ‫ݐ‬, ܶ݉ (° C) is the heat produced by metabolic processes (assumed constant across individuals), and ߝ௜,௧ is the residual term assumed to be distributed as Normal(0, ߪ்ಳ ). The heat transfer coefficient was modeled as ݇௜,௧ = β୩,଴ + ߛ௞,௜ + ߚ௞,ଵ ‫ݔ‬ௐ೎ ,௜ + ߚ௞,ଶ ‫ݔ‬௪௔௥௠,௧ , (Equation 4.4) where β୩,଴ is the intercept; ߛ௞,௜ is an individual random effect assumed to be distributed as Normal(0, ߪఊ ); ߚ௞,ଵ is the coefficient for the effect of centred body weight (‫ݔ‬ௐ೎ ,௜ ); and ߚ௞,ଶ is the effect of a warming trend if ܶ஻೔,೟షభ < ܶ஺ ௜,௧ିଵ(i.e. ‫ݔ‬௪௔௥௠,௧ = 1). A cooling (ܶ஻೔,೟షభ > ܶ஺ ௜,௧ିଵ ) or lack of trend (i.e., ܶ஻೔,೟షభ = ܶ஺ ௜,௧ିଵ , occurring in less than 0.004% of the data) were represented by the intercept when ‫ݔ‬௪௔௥௠,௧ = 0. Body weight was included on the physical assumption that a larger body will take a longer time to heat or cool, and the trend factor was included as cold-adapted species, as habitat specialists, may cool and warm at different rates (Amat-Trigo 2023). 118 I first fit the global model, which included both the body weight (‫ݔ‬ௐ೎ ,௜ ) and warm trend (‫ݔ‬௪௔௥௠,௧ ) effects on k, to determine the best autocorrelation structure (no autocorrelation vs various orders of autoregressive structures) accounting for the strong temporal autocorrelation present in the high frequency (5-second) time series data. Using autocorrelation plots, I determined that a fourth-order autoregressive structure best accounted for the temporal autocorrelation. Then, I fit four candidate models with different combinations of fixed effects ܹ௖ and ‫ ݉ݎܽݓ‬, as well as an intercept only model. For these four models, I used maximum likelihood during preliminary analyses and selected among them using Akaike’s Information Criterion with correction for small sample sizes (AICc; Burnham and Anderson 2002). Models comprising the 95% confidence set for the best model (i.e. whose cumulative AICc weights totaled at least 0.95; Burnham and Anderson 2002) were then refit using maximum restricted log-likelihood. Models were fit using the R package nlme (Pinheiro and Bates 2022), and selected models were visually assessed for fit using residuals vs. fitted plots and QQ-plots. Predicting body temperatures of tagged fish I used these models to generate model-averaged predictions of body temperature from ambient temperature data transmitted by the radio tags. Models including the effect of Wc on k used the body weight (g) of the fish as measured at the time of capture. As the body temperature of the fish was not known at the time of the initial transmission of ambient temperature in a series of detections, initial body temperatures were randomly drawn from a normal distribution with a mean equal to the ambient water temperature at the time of first detection in a series and a standard deviation equal to 1.5°C (this was the maximum standard deviation of temperatures in a study reach at any time during the study). In cases where there were gaps in the detection series lasting less than five minutes, the last estimated body temperature was used to start the next series. In 119 cases where more than five minutes elapsed between detections, another random draw was taken using the new value of the ambient water temperature as the mean of the normal distribution. This process was repeated 100 times for each tagged fish that transmitted enough data for subsequent analyses of behavioural thermoregulation. Estimates of body temperature typically converged in under one hour with a coefficient of variation < 0.05 regardless of the draw for initial body temperature (Supplemental Information SI 4.1). Estimates of body temperature occurring before the coefficient of variation for the 100 simulations dropped below the 0.05 threshold were removed from analyses of behavioural thermoregulation. I note that the heat transfer coefficient k used in the predictive model was parameterized over 5second time steps but, in the radio tagging study, this interval could vary depending on the myriad factors that influence the detection process in radio telemetry under field conditions. While most data points from the selected tags were successfully recorded at the 15-second transmission interval of the radio tags (64% of detections), I did record longer transmission intervals (mostly 30, 45, and 60 seconds [34%], though some longer intervals occurred [2%]). Therefore, my use of the heat transfer equation to convert water temperature measured by the tags to body temperatures assumes that the heat transfer coefficient estimated with data measured at 5-second intervals can be applied to longer time intervals. Heat transfer is a physical phenomenon that is sensitive to interval length and most of the literature investigating this aspect comes from studies in engineering (e.g., Su and Hewitt 2002). To my knowledge, heat transfer across intervals larger than 15-seconds (Pépino et al. 2015) has not been explored in fish telemetry studies. To reduce potential bias in the body temperature predictions resulting from long intervals, intervals lasting more than five minutes were treated as new time series with a new starting value of ܶ஺ (as described above). 120 Behavioural thermoregulation I used three analyses to investigate behavioural thermoregulation by radio-tagged Arctic grayling. All three analyses used metrics of behavioural thermoregulation that were derived from the labbased estimates of ܶௌா் . Because males showed a slightly higher thermal preference than females, ܶௌா் , and thus the resulting metrics, were calculated by sex. Furthermore, all metrics were calculated using each record of ܶ஻ (converted from ambient temperature transmitted by radio tags with equation 3) and associated ambient temperature recorded by temperature data loggers (Supplemental Information SI 4.6). The resulting values were then aggregated as hourly means for each fish for all three analyses. In the first analysis, the absolute deviation between the body temperatures and the thermal preference range, ܾ݀௜,௧ = หܶ஻೔,೟ − ܶௌா் ห, sometimes called thermoregulatory accuracy (with low values representing better thermoregulation; Hertz et al 1993), was regressed against the absolute deviation between ambient water temperatures and the thermal preference range, ݀݁௧ = หܶ஺೟ − ܶௌா் ห, sometimes called environmental quality (with low values representing better quality thermal habitat; Hertz et al 1993): തതതത௜,௧ , ߪఌ ൯, ܾ݀௜,௧ ~Normal൫ܾ݀ ೏್ തതത ܾ݀ത௜,௧ = ߚௗ௕,଴ + ߛௗ௕,௜ + ߚௗ௕,ଵ ݀݁௧ , (Equation 4.5) (Equation 4.6) തതതത௜,௧ is the expected accuracy of thermoregulation for individual ݅ at time ‫ݐ‬, ߪఌ is the where ܾ݀ ೏್ standard deviation of residuals; ߚௗ௕,଴ is the intercept; ߛௗ௕,௜ is a random effect of individual on the intercept which is assumed to be distributed as Normal(0, ߪఊ೏್ ); and ߚௗ௕,ଵ is the slope of the relationship with the quality of the thermal habitat ݀݁௧ . The model also included a third order autocorrelation structure to account for the temporal correlation in the residual errors and was fitted 121 to the data using package nlme (Pinheiro and Bates 2023) in R (R Core Team 2023). The slope of the regression line was compared to the line of thermoconformity (slope = 1) using a t-test. Slopes less than 1 indicate behavioural thermoregulation and the vertical distance between individual observations and the line of thermoconformity is a measure of the degree of thermoregulatory effectiveness of each observation (Blouin-Demers and Nadeau 2005). For the second analysis, I first computed two variables: ‫ݕ‬஻ೄಶ೅ ,௧ , which represented the number of occasions (usually over 15 second intervals; 240 occasions per hour), within hour ‫ ݐ‬that the body temperature of a fish was within ܶௌா் ; and ݊஺ೄಶ೅ ,௧ . This variable represents the number of occasions that ܶௌா் was available in the environment in hour ‫( ݐ‬Blouin-Demers and Weatherhead 2001; Blouin-Demers and Nadeau 2005). I then used a binomial generalized linear mixed model to assess whether the proportion of time a fish spent within ܶௌா் – termed thermal habitat exploitation and represented by ‫ܧ‬௑ – varied as a function of the quality of the thermal habitat (݀݁௧ ): ‫ݕ‬஻ೄಶ೅ ,௧ ~Binomial(݊஺ೄಶ೅,೟ , ‫ܧ‬௑,௧ ), (Equation 4.7) Logit൫‫ܧ‬௑,௧ ൯ = ߚ௬,଴ + ߛ௬,௜ + ߚ௬,ଵ ݀݁௧ , (Equation 4.8) where ‫ܧ‬௑,௧ is the proportion of time in hour ‫ ݐ‬in which an the body temperature of an individual remained within their thermal preference range ܶௌா் on occasions that it was available in their environment; ߚ௬,଴ is the intercept; ߛ௬,௜ is a random effect of individual on the intercept which is assumed to be distributed as Normal(0, ߪఊ೤ ); and ߚ௬,ଵ is the slope of the effect of ݀݁௧ on logit൫‫ܧ‬௑,௧ ൯. The model was fitted to the data using penalized quasi-likelihood methods in the package MASS (Venables and Ripley 2002). The third analysis focused on the index of behavioural thermoregulatory effectiveness ‫୼ܧ‬ (Christian and Weavers 1996; Blouin-Demers and Weatherhead 2001; Blouin-Demers and Nadeau 122 2005). This index is computed as the difference between ݀݁ at time ‫ ݐ‬and the ܾ݀ by individual ݅ at time ‫ݐ‬: ‫୼ܧ‬,௜,௧ = ݀݁௧ − ܾ݀௜,௧ (Equation 4.9) Index ‫୼ܧ‬,௜,௧ quantifies thermoregulatory effectiveness as either positive (higher values indicating increasing thermoregulatory effectiveness as it means the difference between the body temperature and preferred range is smaller than the difference between the ambient temperature and preferred range) or negative (lower values indicating increasing avoidance of preferred temperatures) departures from zero (thermoconformity). This index is a rescaled version of the original formulation of index E for effective behavioural thermoregulation (Hertz et al. 1993; BlouinDemers and Weatherhead 2001). Index ‫୼ܧ‬,௜,௧ was modeled as a function of a suite of covariates using generalized additive mixed models (GAMMs) with the form ‫୼ܧ‬,௜,௧ ~Normal൫‫ܧ‬ത୼,௜,௧ , ߪఌಶ ൯, (Equation 4.10) ‫ܧ‬ത୼,௜,௧ = ߚா,଴ + ߛா,௜ + ݂൫‫ݔ‬௔௖௧೔,೟ ൯ + ݂൫‫ݔ‬௣௔௧௖௛೟ ൯ + ݂൫‫ݔ‬௛௢௨௥೟ ൯ × ‫݊݋ݏܽ݁ݏ‬௧ + ߚா,ଵ ‫ݔ‬௕௢ௗ௬௖௢௡೔ + ߚா,ଶ ‫ݔ‬௦௘௔௦௢௡೟ , (Equation 4.11) where ‫ܧ‬ത୼,௜,௧ is the expected effectiveness of thermoregulation for individual ݅ at time ‫ݐ‬, ߪఌಶ is the standard deviation of residuals; ߚா,଴ is the intercept, ߛா,௜ is the random effect of individual on the intercept assumed be distributed as Normal (0, ߪఊಶ ); ݂(. ) denotes smoothing functions of hourly mean activity (ܽܿ‫ )ݐ‬exhibited by individual ݅ at time ‫ݐ‬, the patchiness (‫ܿݐܽ݌‬ℎ) of the thermal environment and hour of the day (ℎ‫ )ݎݑ݋‬at time ‫ ;ݐ‬and ߚா,ଵ represents the effect of the body condition (ܾ‫ )݊݋ܿݕ݀݋‬of individual ݅. While conducting preliminary data visualization, I noticed a change in the hourly patterns of ‫ ୼ܧ‬after September 4, so I included a categorical variable ‫݊݋ݏܽ݁ݏ‬௧ 123 as an interaction on ℎ‫ݎݑ݋‬௧ , which and indicates whether an observation at time ‫ ݐ‬occurred before and after a shift to cooler temperatures on September 4. A total of 32 candidate GAMMs with different combinations of covariates were fit and selected by Akaike’s Information Criterion with correction for small sample sizes (AICc; Burnham and Anderson 2002). The models included in the 95% confidence set for the best model were used to create model-averaged predictions of covariate effects on thermoregulatory effectiveness. GAMMs were fit using 20 knots for each smoothed term with the R package mgcv (Wood 2011). The thin plate regression spline was used for all smoothed terms except for hour, which used the cyclic cubic regression spline (Wood 2011). Model averaging was done using the package MuMIn (Barton 2023). To account for strong temporal autocorrelation in the model, GAMMs were fit using a third-order autoregressive structure (Pinheiro and Bates 2023) and assessed for fit using residual vs fitted plots and QQ-plots. Results Thermal preference and heat transfer The thermal preference experiments revealed a median thermal preference (ܶ௉ோாி ) of 11.9 °C (± 4.1 °C median absolute deviation [MAD]; n = 10) with a range (ܶௌா் ) between 10.1 °C and 13.0 °C (Table 4.2). Median ܶ௉ோாி of females was lower (10.2 ± 2.7 °C MAD; n = 6) than in males (13.5 ± 4.3 °C MAD; n = 4), but the difference was not statistically significantly (t = 1.1 p = 0.3; DF = 4.3). The heat transfer model that had the most support included both body weight and the heating/cooling trend (‫ )݉ݎܽݓ‬effects on the heat transfer coefficient ݇ (Table 4.3). A second model 124 including only ‫ ݉ݎܽݓ‬was also included in the 95% confidence set (Table 4.3). Predictions from the best model when body and water temperature differed by 3°C in a hypothetical scenario revealed that it would take between 10-15 min for the body and water temperature to equilibrate, with smaller individuals equilibrating more quickly than larger ones (Figure 4.2). The rates of heating and cooling were similar, with slightly higher efficiency when cooling (Supplemental Information SI 4.4). Due to metabolic heat production quantified by ܶ௠ , slightly higher body temperature after equilibration were achieved in predictions for larger fish than for small fish (Figure 4.2). Behavioural thermoregulation തതതത௜,௧ = ߚௗ௕,଴ + ߛௗ௕,௜ + ߚௗ௕,ଵ ݀݁௧ ; Hertz et al. 1993), revealed that The estimated slope for Eq. 4.6 (ܾ݀ the index quantifying body temperature deviations (ܾ݀) by Arctic grayling increased (i.e. body temperature deviated further from ܶௌா் ) as the index quantifying the habitat temperature deviations (݀݁) increased (i.e. the thermal habitat deviated further form ܶௌா் ; Figure 4.3). However, the slope describing the change in ܾ݀ per unit change in ݀݁ was only 0.473, which is significantly smaller than the slope of 1 associated with thermoconformity (t = 22.9; df = 3523; P < 0.001). This indicates that Arctic grayling were able to behaviourally maintain body temperatures closer to ܶௌா் as habitat temperature deviations increased more than would be expected if they did not exhibit behavioural thermoregulation. However, as habitat temperature deviations increased, exploitation of preferred thermal habitats (index ‫ܧ‬௑ ), when they were available, decreased (Figure 4.4; slope = -0.879; t = -9.43; ݂݀ = 1590; P < 0.001). Plots of ܶ஻ , ܾ݀, ‫ ୼ܧ‬, and ‫ܧ‬௑ across all study hours temporal patterns of Arctic grayling thermoregulatory behaviours, with the most pronounced patterns visible along the diel cycle (e.g. Figure 4.5; all plots are available in Appendix B). Across fish with sufficiently long time series, a 125 pattern can be observed where effective thermoregulatory behaviours (positive ‫ ) ୼ܧ‬more commonly occurred during the day in the warm part of the season. However, after approximately hour 850 (September 4), when temperatures dropped and were consistently below ܶௌா் for most fish, effective thermoregulation became more common at night. Heterogeneity of thermal habitats in space negatively covaried with both diel and seasonal temporal heterogeneity, though correlation was low (-0.02 with diel and -0.07 with seasonal). Thermal exploitation (Index ‫ܧ‬௑ ) was variable during the warm season, but for fish with sufficiently long time series, index ‫ܧ‬௑ dropped to zero for most observations after September 4. The top ranked thermoregulatory effectiveness (index ‫ ) ୼ܧ‬GAMM included the interaction between ℎ‫ ݎݑ݋‬and ‫ ݊݋ݏܽ݁ݏ‬and held 42% of the AICc weight (Table 4.4). These terms were included in all seven models in the 95% confidence set, but there was substantial uncertainty in the other terms which were not included across all top models. Condition factor (ܿ‫ )݊݋‬was present in the second-ranked model and two other models of the 95% confidence set, whereas patchiness (‫ܿݐܽ݌‬ℎ) and activity (ܽܿ‫ )ݐ‬were included in only two models (Table 4.4). Model-averaged predictions revealed that thermoregulatory effectiveness (index ‫ ) ୼ܧ‬was higher when ambient temperatures were cooler in the fall (after September 4; Figure 4.6). Hour of day showed a strong nonlinear effect with thermoregulatory effectiveness, with fish being more effective thermoregulators in the mornings and afternoons during the warm season (before September 4th), and strong thermoregulators in the mornings during the cool season (after September 4th; Figure 4.6A). The highest point of thermoregulatory effectiveness occurred at approximately 14:45 in the warm season, which corresponded with the warmest temperatures recorded over the diel cycle, and this reversed in the cool season, where thermoregulatory effectiveness was lowest during the warmest part of the day. Thermoregulatory effectiveness index 126 ‫ ୼ܧ‬showed a weakly negative linear relationship with condition factor, where fish with higher condition factors showed a marginally lower thermoregulatory effectiveness (Figure 4.6B). The patchiness of the thermal habitat also showed a weak nonlinear effect on thermoregulatory effectiveness, with fish being marginally less effective thermoregulators in intermediately patchy thermal environments (Figure 4.6C). Index ‫ ୼ܧ‬was weakly linearly related to mean hourly activity, with thermoregulation being slightly more effective at low values of activity and decreasing as activity increased (Figure 4.6D). Discussion In this chapter, I quantified the thermal preference of Arctic grayling and demonstrated they use behaviour to maintain their body temperatures within preferred temperatures during their summer feeding window. I found that despite the reduction in thermal exploitation as thermal habitats moved further away from ܶௌா் , the index of effective behavioural thermoregulation increased. In other words, Arctic grayling were careful to maintain their body temperatures close to ܶௌா் even in poor-quality thermal habitats. I tested the hypotheses that thermoregulatory effectiveness in Arctic grayling was driven by diel hour, condition factor, and the patchiness of the thermal habitats, and that these strategies were energetically costly (i.e. required more activity) to maintain. I found that the effectiveness of thermoregulation was related more strongly to some of the effects tested than others. In particular, I found that the heterogeneity of thermal habitats over time (hour of day and seasonal change) was a much stronger predictor of thermoregulatory effectiveness than heterogeneity in space, and that that condition factor and activity had minor effects on effectiveness (though the range of condition factor [0.069 – 0.209] may have been limiting compared to other measures). 127 The median ܶ௉ோாி of 11.9 °C and ܶௌா் of 10.1 °C - 13.0 °C revealed by the shuttlebox experiments for the Anzac River Arctic grayling presented a reasonable thermal preference range for montane fluvial fishes in an Arctic watershed. This range was lower than the well-defined ܶௌா் for more temperate latitudes salmonids (e.g., ܶ௉ோாி has been identified as 14.9 – 15.0 °C for westslope cutthroat trout [Oncorhynchus c. lewiski], 16.0 °C for brown trout [Salmo trutta], and 14.8 °C for rainbow trout [Oncorhynchus mykiss]; McMahon et al. 2008; Larsson 2005; Macnaughton et al. 2021), was similar to that of other Arctic salmonids (e.g. two populations of Arctic char [Salvelinus alpinus] were estimated as having ܶ௉ோாி at 10.8 °C and 11.8 °C; Larsson 2005), and fell within estimates of habitat parameters from other Arctic grayling populations. Habitat suitability indices estimate optimal thermal habitats of 6 – 16 °C for adult Arctic grayling (Hubert et al. 1985; Stewart et al. 2007) and 5 – 14 °C across all life stages (Larocque et al. 2014). A site occupancy study in the Parsnip River watershed found a high probability of occurrence at temperatures of 10.9 °C (with > 0.75 probability of occurrence within 8.7 – 14.2 °C; O’Connor 2023). My estimate of ܶௌா் also agrees with the findings of watershed-scale habitat use occurring at 11.3 °C (11.1 – 17.1 °C) found in chapter 2 of this dissertation and the reach-scale estimates of 12.5 °C (9.2 – 16.0 °C) from chapter 3. While one individual (the largest male) used in the shuttlebox experiments showed an unusually high ܶௌா் of 17.6 – 18.8 °C, the mean ܶௌா் across all experiments also supported the avoidance of stressful temperatures above 17.2 °C (Wojcik 1955; Cahill 2015; Carillo-Longoria 2023) and was well below estimates of incipient lethal temperatures of 20 – 25 °C (Wojcik 1955; Lohr et al. 1996; Cahill 2015; Tingley et al. 2022). The estimated ܶ௉ோாி in this study was lower than an experimental study on Arctic grayling juveniles which found a ܶ௉ோாி of 16.7 °C (Hawkshaw 2011). The ten fish which were successfully trialed for thermal preference showed 128 TPREF ranging from 7.2 – 18.4 °C, which demonstrates that TPREF is a variable metric within populations subject to the same availability of thermal habitats. Population-level variability in thermal responses can be challenging for conservation planners dealing with species at the southern peripheries of their distribution (Dressler et al. 2023). While thermal preference is known to be a relatively fluid metric that can change both seasonally and within the lifespan of an individual as life history goals change (Bloomfield et al. 2022), studies on other Arctic fish species suggests that ܶ௉ோாி may be relatively stable throughout the open-water season (Norwegian Arctic char had a stable ܶ௉ோாி from 11.5 - 11.8 °C from spring to fall; Mortensen et al. 2007). As an Arctic species at the southern periphery of their distribution in B.C., it is not known if or at what rate and variability the ܶ௉ோாி changes seasonally in my study population (though they likely co-vary). Shuttlebox experiments occurred from August 25 – September 10, 2022, approximately 24% of the duration of the in-situ radio tagging study, and acclimation temperatures (mean river temperatures for 2 weeks prior to the beginning of each shuttlebox experiment) were 13.3 °C (±1.1 °C; Supplemental Information SI 4.3). As such, it is likely that the ܶ௉ோாி did not change significantly at the individual level over this period. My estimates of the heat transfer coefficient k were a function of both body weight (on the physical assumption that larger bodies are slower to heat or cool) and the direction of the temperature difference between body temperatures and ambient water temperatures (on the assumption that cold-water adapted species may warm and cool at different rates) over each time interval. Body weight of fish was found to be important to estimating heat transfer in Mozambique tilapia (Tilapia mossambia; Stevens and Fry 1974), bluegill (Lepomis macrochirus), blue tilapia (Oreochromis aureus), Nile tilapia (Oreochromis niloticus; Fechhelm and Neill 1982), largemouth bass (Micropterus salmoides; Weller et al. 1984), and brook char (Salvelinus fontinalis; Pépino et al. 129 2015). While Pépino et al. (2015) found that the magnitude of the temperature differential was important for brook char (a factor which is accounted for in the heat transfer equation), they found that the direction of the heating or cooling trend was not important to the rate of heat transfer. Weller et al. (1984) found a significant effect of body weight on whether largemouth bass were heating or cooling, though this is a warm-water species (thermal preferences between 28 – 29°C in juveniles with higher estimates in adults; Diaz et al. 2007) and largemouth bass are among a subset of species known to thermoregulate in both directions (Amat-Trigo et al. 2023). While the trend variable (warm) was found to be important to explaining the heat transfer coefficient k, there was only a small difference between k under heating and cooling conditions. However, the higher rate of heat transfer found when cooling, even slight, was an interesting finding given Arctic grayling are a cold-water adapted species and is worth exploring further. The regression analysis of the accuracy of thermoregulation, ܾ݀, against the thermal quality of the habitat, ݀݁, showed that Arctic grayling thermoregulate more effectively as their environment moves further away from their ܶௌா் . However, the low ܶௌா் exploitation at only moderately increased values of ݀݁ suggest that this behaviour is mediated by other tradeoffs presented by their environment. It is well understood in habitat ecology that an animal will exploit trade-off habitats if the benefits of doing so outweigh the costs (Veech 2021), a behaviour that has also been observed in thermoregulating freshwater fishes (McCollough et al. 2009; White et al. 2019; Amat-Trigo 2023). In this study population, it is likely that feeding opportunities were drawing Arctic grayling away from their preferred thermal habitats during the day or Arctic grayling were occupying warmer waters temporarily as a means to aid the efficiency of digestion during summer rearing (Armstrong et al. 2021). This effect can be observed in bioenergetics models for bull trout (Salvelinus confluentus) that have shown that optimal energetic intake occurs at 16 °C (Mesa et al. 130 2013), which is higher than estimates of their thermal preference (~ 12 °C; Stewart et al. 2007). This behaviour likely applies to Arctic grayling which are known to distribute themselves relative to optimal energy intake during the summer months (Hughes 1992a,b). As Arctic grayling are highly reliant on vision to feed (McPhail 2007; Stamford et al. 2017), these behavioural tradeoffs are only likely to occur during the daylight hours when their preferred prey of terrestrial drift can be spotted. Indeed, a clear diel pattern was found in the analysis of thermoregulatory effectiveness ‫ ୼ܧ‬. In the summer (i.e. the warm season before September 4), thermoregulatory effectiveness was at its highest point as afternoon temperatures peaked at approximately 2:45 PM. This coincided with approximate solar noon in the Anzac River valley (the approximate flight times of the drone surveys in Table 3.1 of Chapter 3). Thermoregulation decreased from this peak in the mornings and early evening hours, times during which both visual feeding was possible and corresponded with insect hatches. These times likely represent the tradeoffs in foraging and thermoregulating behaviours. A second, smaller peak in thermoregulatory effectiveness in the morning hours (approximately 9:00 AM) suggests that Arctic grayling can effectively navigate this trade-off during the cool morning hours. Interestingly, overnight hours showed a pattern of negative ‫ ୼ܧ‬, indicating avoidance of preferred thermal habitats (BlouinDemers and Nadeau 2005). It is possible that Arctic grayling not experiencing diel temperature stress simply thermoconform at night rather than seeking out warmer temperatures. This behaviour may also be a result of nighttime predator avoidance (in that seeking out preferred thermal habitats would expose them to predation risk) or to aid in digestion when not feeding (Brett 1970). Thermoconformity at night further suggests that Arctic grayling may be single-direction thermoregulators that only invest energy into cooling, but may be comfortable thermoconforming at sub-ܶ௉ோாி water temperatures. Other cold- and cool-water salmonid species have demonstrated 131 single-direction thermoregulation during the summer months (e.g. bull trout Salvelinus confluentus, brook trout Salvelinus fontinalis, sockeye salmon Oncorhynchus nerka, and cutthroat trout Oncorhynchus clarkii [Donaldson et al. 2009; Howell et al. 2012; Hodge et al. 2017; Goyer et al. 2014; Hitt et al. 2016; Amat-Trigo et al. 2023]). As Arctic grayling are cold-water specialists, this finding seems reasonable and further research that explicitly measures the directions of temperature deviations (as opposed to absolute deviations used in this study) would be useful. With the start of the cool season after September 4, morning thermoregulation became more pronounced, and there was not a second peak in thermoregulatory effectiveness in the afternoon (though overall thermoregulatory effectiveness was higher in this season). Across the sampling season, water temperatures were highest in the study reaches around 15:00, the same time that a reverse trend was seen in thermoregulatory effectiveness (i.e. the time of highest effectiveness in the warm season was the lowest in the cold season; Figure 4.6). This may be attributed to differential responses to unfavorably warm and unfavorably cool temperatures in cold-water adapted Arctic grayling. While unfavorably cold waters can suppress metabolic activity, extreme heat can quickly become lethal, and even the coldest temperatures recorded in the study reaches were within habitat suitability estimates for Arctic grayling (Larocque et al. 2014). Studies which use the E indices of behavioural thermoregulation (Hertz 1993; Blouin-Demers and Weatherhead 2001; Blouin-Demers and Nadeau 2005) have been predominantly focused on reptile populations and to my knowledge only one other study (being conducted within our working group) has quantified thermoregulatory effectiveness in free-ranging fishes (Dextrase 2024). It is possible that the thermoregulatory effectiveness results are a function of environmental structure; favourable thermal habitats may be available during certain hours of the diel cycle in unsuitably shallow or high-flow waters. As thermal habitats in rivers follow a strong hourly and daily cycle 132 (Caissie 2006), all available habitats will pass through the suitable thermal range at different times of the day. Available river temperatures were defined using the FLIR rasters created in Chapter 3 and these were cut to the wetted portion of the images. While anecdotal evidence suggests that Arctic grayling will at times use shallow habitats near the river’s edge (J. Richert, West Moberly First Nations, pers. comm.), it is not likely that they were regularly used for thermoregulatory purposes as the monitored fish tended to remain in deeper pool habitats most of the time. High variation in index ‫ܧ‬௑ when ݀݁ = 0 (when most of the data observations were available) indicate that Arctic grayling were free to navigate foraging tradeoffs in high-quality thermal habitats. However, insights from index ‫ܧ‬௑ are limited as it does not take into account thermoregulatory behaviour during periods where the environment is the most thermally challenging (i.e., when ݀݁ > 1; Blouin-Demers and Nadeau 2005). I note that these results may have been influenced over fine-scales by the low temperature resolution of 0.4 °C of the transmitters, which may have obscured inferences into thermoregulatory effectiveness when ܶௌா் and ݀݁ were close together. Contrary to my hypothesis that behavioural thermoregulation is energetically costly (i.e. requires fish to be more active) for Arctic grayling, thermoregulatory effectiveness did not continually increase with higher activity, but rather were slightly higher at low values of activity and decreased modestly as activity continued to increase. The decrease in thermoregulatory effectiveness at high values of activity could have been a result of a few factors. It is possible that small energetic expenditures do contribute to thermoregulatory effectiveness up to an extent after which expending more energy simply does not provide a significant thermoregulatory return. While the energetics of thermoregulation have not been well-explored in fishes despite the advances in activity-sensing radiotelemetry (Lennox et al. 2023), they would be theoretically expected to be costly to maintain 133 (Huey and Slatkin 1976; Blouin-Demers and Nadeau 2005; Parlin and Schaeffer 2022). However, a study in lizard populations found that while increased energy did not provide a thermoregulatory return, they experienced a slower growth rate (Brewster et al. 2013), suggesting that there is an indirect diminishing return on energy expenditures. This finding could also be a signal derived from the feeding tradeoffs described above; availability of forage for Arctic grayling in fluvial environments is a function of flows (Hughes and Dill 1990; Hughes 1992a; Hughes 1992b; Hughes 2000) and Arctic grayling must expend more energy to maintain their position when feeding in high flow environments relative to areas of flow refugia. It is also possible that high values of activity were associated with other costly behaviours, such as avoiding predators or moving across shallow riffles between feeding pools and, in either case, thermoregulation would not be the immediate goal. In the latter case, movements between habitats could certainly have the end goal of occupying more favourable thermal habitats, though it is not known how far Arctic grayling might distribute within their reach-scale habitats in search of suitable temperatures. Juveniles of other salmonid species have been documented using horizontal diel migration in thermoregulation on the scales of 300 – 1700 m (Armstrong et al. 2015). The unexpected results about energy expenditures could also be explained by the how the tags measured3 fish activity. The majority of activity detections in this study were from small movements (95.5% of jerk detections were between 0 – 0.2 on a scale to 1.5) and higher measurements of activity were limited in the data. During the snorkeling study in chapter 3, Arctic grayling that were observed feeding would hold their position in the current and quickly dart up to forage before returning to their original positions. These burst movements happened in a matter of seconds; events lasted a relatively small proportion of the radio tag’s 15-second averaging cycle for activity and recorded values may have been smoothed by the remainder of the period’s 134 relatively small movements. It is possible to interpret the activity values recorded near zero as having some motion along the sway axis while very little activity would register along the surge or heave axes during periods when Arctic grayling were holding their positions in the current (sensu Qasem et al. 2013). The placement of the tag may have contributed to this as well; on smallbodied fish such as Arctic grayling it had to be placed in the middle of the body under the dorsal fin, but an ideal placement for measuring activity would have been on the tail where the most sway would be detected. It would be possible to obtain higher-resolution activity data if the tags were retrieved and manually downloaded, but over the course of this study over 70 Arctic grayling were sampled without a single recapture and the opportunity was not presented. The effect of the patchiness of thermal habitats was an uncertain and weak contributor to thermoregulatory effectiveness. Much of the theory underpinning the importance of thermally heterogeneous environments to behavioural thermoregulation was developed using studies on terrestrial ectotherms (Huey and Slatkin 1976; Blouin-Demers and Weatherhead 2001; BlouinDemers and Nadeau 2005; Sears and Angiletta 2015; Sears et al. 2016), and studies on thermal patchiness in rivers have been focused on identifying specific areas associated with cold-water refugia (Casas-Mulet 2020; Kuhn et al. 2021). To-date no studies have attempted to explore these dimensions in fish populations. As the metrics of thermal availability within the reaches were defined using both the FLIR rasters and temperature loggers used in chapter 3, there is some uncertainty as to whether all of the heterogeneity within the thermal habitats of the study reaches was captured using my methods. For example, in the 1 – 2 km study reaches, patchiness was calculated over the total reach; whether fish at the head of the reach would consider temperatures in the middle of the reach was not known. While the snorkeling study in chapter 3 found some redistribution of abundance between pool habitats over 2-week intervals, this dynamic was not 135 well-understood at the hourly scale over which the effectiveness of thermoregulation was modeled. As the radio telemetry setup used in this study was programmed for biologging rather than positioning, the exact locations of the radio tags through time were not available to explore this dimension further and this may have been the source of some uncertainty in the models. While the theory that heterogeneity in the thermal environment is a key consideration for effective behavioural thermoregulation is likely still relevant to aquatic systems, my measurements of thermal patchiness may have been inadequate to capture this across appropriate scales. Alternatively, the strong diel effects in rivers and the high specific heat of water may have masked this variable in my models, or thermoregulation in fluvial Arctic grayling may simply be more driven by temporal heterogeneity than spatial heterogeneity. My FLIR rasters were captured following existing recommendations for thermal drone studies in rivers and surveys were conducted at solar noon at which the thermal heterogeneity of river habitats is most visible by thermal sensors (Dugdale 2016). Rasters of temperature distributions were interpolated based on instream temperature loggers and assume that the spatial heterogeneity of the thermal habitats were relatively static and a function of flowing water through river morphology. More frequent drone flights across more varied times of day, and in particular night surveys, may have captured further heterogeneity in stream temperature distributions to allow for a finer-scale examination of the spatial thermal heterogeneity in the study reaches (though shadows and other factors would have influenced these as well). Examining the relationship and magnitude between spatial and temporal heterogeneity thermal in rivers would be a useful area for future research which could further clarify the relative importance of spatial and temporal heterogeneity. The final metric explored in the analysis, condition factor, showed a weakly negative effect on the effectiveness of behavioural thermoregulation. This is likely an effect that can be explained using 136 the same physical rationale as modeling coefficient k on body weight; fish with a higher condition factor are likely to have more thermal inertia (Andrare et al. 2015) and may not have to regulate their body temperatures as carefully as fish who are more susceptible to the magnitude and direction of changes in their thermal environment. Because body condition is an integrated metric that considers both length and weight, it is possible that fish with low condition factors may have to be more careful at thermoregulating as they might not have as much energy to spare as fish with higher condition factors (though this also interplays with other physiological states; Brosset et al. 2023). Together, these results depict a complex daily interplay among the tradeoffs of feeding, spatiotemporal thermal heterogeneity, avoiding thermal extremes and apparent thermoconformity to cold temperatures, condition factor, and digestion and growth. As Arctic grayling navigate their diel thermal experience during their summer feeding window, their responses to temperature extremes and stressors which exceed their thermal preference range will be proportional to both the magnitude of duration of thermal stressors, and the spatiotemporal heterogeneity of adequate recovery temperatures (Farrell et al. 2008; Rezende et al. 2014; Wolkovich et al. 2014). As climate change forces more frequent and more severe extreme temperature events (like the heat dome event across British Columbia in the summer of 2021), Arctic grayling use of behavioural thermoregulation to navigate the tradeoffs presented by their environment will be imperative to their long-term survival. Thermoregulatory behaviours can help buffer extreme temperatures over the short term and will be an important component in climate change adaptation, but evidence shows that these behaviours can ultimately weaken selective pressures on populations, slowing population-scale adaptation to climate change over generations (Buckley et al. 2015). Further, stress responses over mild temperature extremes can improve heat tolerance while severe extremes 137 can diminish it (Rodgers and Isaza 2022). In light of an uncertain future, habitats which support the thermal preference range of Arctic grayling and continued exploration of how individuals exploit them over space and time, will continue to become more important to conservation planning in the Parsnip River watershed and beyond. 138 Chapter 4 Tables Table 4.1. Experimental design for the shuttlebox heat transfer experiments. Each of the five experimental trials (shifts) exposed an individual Arctic grayling to a warm (+) or cool (-) temperature differential. The initial direction of the temperature shift pattern was determined by a randomized start in either the cool tank or the warm tank. As trials were conducted, the magnitude of the temperature differential gradually decreased from a starting differential of 2.0 °C to a final differential of 1.0 °C. Starting in cool tank Starting in warm tank Shift Direction °C Shift Direction °C 1 2 3 + - 2.0 2.0 1.5 1 2 3 + + 2.0 2.0 1.5 4 + 1.0 4 - 1.0 5 - 1.0 5 + 1.0 139 2022-08-28 2022-08-29 2022-08-31 2022-09-02 2022-09-03 2022-09-05 2022-09-06 2022-09-08 2022-09-09 2022-09-10 2 3 4 5 6 7 8 9 10 Date 1 Trial 19:36 16:30 15:01 19:07 07:50 18:10 15:15 08:35 20:18 09:35 Capture time 10.2 10.2 10.6 10.7 9.1 13.9 14.4 9.1 12.2 9.8 Capture temp (°C) 62 61 60 59 58 57 56 55 54 53 Fish no. M F F M M F F M F F Sex 355 304 345 330 362 293 343 353 332 312 Len. (mm) TPREF ± TSET : 500 450 500 425 540 350 500 660 490 360 Wt. (g) 10.1 10.3 8.0 9.2 13.1 17.6 6.2 7.4 7.0 10.5 11.6 q25 11.9 12.6 9.0 13.9 14.5 18.4 7.2 9.0 7.6 11.3 12.6 q50 TOCC (°C) 13.0 15.2 9.9 17.0 15.5 18.8 9.6 9.7 8.4 12.8 13.2 q75 140 Table 4.2. Summary table for the 10 Arctic grayling used in the shuttlebox thermal preference experiments. The thermal preference (ܶ௉ோாி ) is given as the median (q50) of occupied temperatures (ܶை஼஼ ) during the twenty hours of data collection and the thermal preference range (ܶௌா் ) is given as the temperatures in the 25 th (q25) and 75th (q75) quantiles. Table 4.3. AICc statistics for the non-linear mixed effects models of the heat transfer coefficient ݇ against centred body weight (Wc) and whether the body temperature of the fish was warming or cooling relative to ambient water temperatures (trend [warm]). Presented are the log likelihood, the number of parameters in the model (K), the AICc score, the difference in AICc between each model from the top model (ΔAICc), and the respective model weights (AICcwt), and the 95% confidence set (highlighted in bold). A random effect was included in the model for each fish used in the experiment. Fixed effects model k ~ Wc + trend (warm) k ~ trend (warm) k ~ Wc k ~ no effects K AICc logLik ΔAICc AICcwt 6 5 5 4 -15588.89 -15583.96 -15583.00 -15578.24 7800.44 7796.98 7796.50 7793.12 0.00 4.93 5.89 10.65 0.88 0.07 0.05 0.00 141 11 11 12 12 -1730 13 -1728.9 14 -1728.8 8 -1746.8 9 -1746.4 10 -1745.5 11 -1745 10 -1746.2 11 -1745.8 12 -1745.3 13 -1744.4 7 -1756.1 8 -1756 9 -1755.2 9 -1755.5 10 -1755.1 EΔ ~ hour x season + patch EΔ ~ hour x season + activity EΔ ~ hour x season + con + patch EΔ ~ hour x season + con + act EΔ ~ hour x season + patch + act EΔ ~ hour x season + con + patch + act EΔ ~ hour EΔ ~ hour + con EΔ ~ hour + patch EΔ ~ hour + con + patch EΔ ~ hour + act EΔ ~ hour + con + act EΔ ~ hour + patch + act EΔ ~ hour + con + patch + act EΔ ~ season EΔ ~ season + con EΔ ~ season + patch EΔ ~ season + act EΔ ~ season + con + patch -1729 -1730 -1729 -1730 -1730 9 10 EΔ ~ hour x season EΔ ~ hour x season + con logLik K Model 3528.367 3529.135 3530.293 3514.933 3526.258 3527.958 3512.43 3513.718 3514.693 3510.894 3510.972 3512.081 3483.55 3483.869 3485.628 3509.603 3481.00 3481.75 3482.76 3478.91 3480.71 AICc 49.46 50.22 51.38 36.02 47.35 49.05 33.52 34.81 35.78 31.98 32.06 33.17 4.64 4.96 6.72 30.69 2.09 2.84 3.85 0.00 1.80 ΔAICc 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.04 0.04 0.01 0.00 0.15 0.10 0.06 0.42 0.17 AICcwt 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.95 0.99 1.00 1.00 0.74 0.85 0.91 0.42 0.60 cAICcwt 142 Table 4.4. AICc statistics for the candidate GAMMs modeling the index of thermoregulatory effectiveness ‫ܧ‬௱ against hour of the day (hour), condition factor (con), activity (act), season (before and after September 4), and patchiness (patch) of the thermal habitat. Presented are the log likelihood (logLik), the number of parameters in the model (K), the AICc score for small sample sizes, the difference in AICc between each model from the top model (ΔAICc), the respective model weights (AICcwt), and the cumulative model weights (cAICcwt) for selecting the 95% confidence set for the best model (highlighted in bold) 6 7 8 9 8 9 10 -1773.9 11 -1773.3 EΔ ~ 1 EΔ ~ con EΔ ~ patch EΔ ~ con + patch EΔ ~ act EΔ ~ con + act EΔ ~ patch + act EΔ ~ con + patch + act -1773.9 -1775 -1774.6 -1775.6 -1775.1 -1774.5 10 -1755.4 11 -1754.6 12 -1754.6 EΔ ~ season + con + act EΔ ~ season + patch + act EΔ ~ season + con + patch + act 3568.736 3567.855 3565.828 3566.093 3567.158 3563.217 3564.287 3564.941 3530.834 3531.276 3533.204 89.82 88.94 86.92 87.18 88.25 84.31 85.38 86.03 51.92 52.36 54.29 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 143 Figure 4.1. Distribution of temperature loggers and radio receivers used in this study. Chapter 4 Figures 144 145 Figure 4.2. Predictions of heat transfer in Arctic grayling along a body weight gradient as exposed to a heating (12 – 15 °C) and cooling (15 – 12 °C) temperature shift. Smaller fish were quicker to heat or cool than larger fish, but larger fish equilibrated at a slightly higher temperatures due to higher metabolic heat production. The rate at which body temperatures equilibrated was similar, with modeling indicating slightly more efficient cooling than heating. Figure 4.3. Regression analysis of the index of accuracy of thermoregulation, ܾ݀, against the index of thermal quality of the habitat, ݀݁. The dashed line represents the line of thermoconformity (slope = 1). Points above this line indicate ineffective thermoregulation and points below this line indicate effective thermoregulation. The regression line of the points falls below the line of thermoconformity and indicates thermoregulatory behaviours. Plot points are displayed with transparency; darker points indicate more overlap in the data. 146 147 Figure 4.4. Regression analysis of the degree of hourly thermal exploitation ‫ܧ‬௫ , or the proportion of time Arctic grayling occupied their thermal preference range (ܶௌா் ) when it was available in their environments. 148 Figure 4.5. An example plot from Arctic grayling no. 15 showing (A) body temperature in blue against median water temperatures in red with the black horizontal lines denoting ܶ௉ோாி (solid) and the upper and lower bounds of ܶௌா் (dashed) against study hour, (B) the accuracy of body temperatures with respect to ܶௌா் ܾ݀, (C) index ‫ܧ‬௱ as positive or negative departures from zero (thermoconformity) with positive values representing thermoregulatory behaviours and negative values denoting avoidance of preferred temperatures and day or night period by color (a general trend from day to night thermoregulation can be seen starting around hour 850 on September 4), and (D) the index of thermal exploitation (‫ܧ‬௑ ) against study hour by day and night. Plots for all fish with sufficient data can be found in Appendix B. Figure 4.6. Model-averaged prediction plots for the GAMM analysis of the thermoregulatory effectiveness (index ‫ ) ୼ܧ‬against (A) mean hourly activity, (B) body condition, (C) the patchiness of the thermal habitat, and (D) hour of the day by season. The dotted horizontal line represents the line of thermoconformity, and the shaded areas bounded by the dashed lines represent the standard error of the predictions. 149 Chapter 5 - Discussion Freshwater biodiversity is declining at a rate disproportionate to other ecosystems, suffering an 83% loss in populations from 1970 to 2014 (Dudgeon et al. 2006; Reid et al. 2018). These environments make up only 2.3% of the world's water supply (Reid et al. 2018), most of which is held in lakes and reservoirs and only a small portion of which flows through rivers. Despite being in the minority of aquatic environments worldwide, rivers serve critical ecosystem functions. Rivers transport water and nutrients, serve as migratory corridors for both fish and wildlife (Sutherland et al. 2015), can be seasonally important habitats to both fluvial and adfluvial species (Hagen et al. 2020), provide flood control function (Anderson et al. 1996), and are home to resident fluvial species. These highly dynamic ecosystems are characterized by shifting abiotic processes along daily, seasonal, annual, and interannual temporal scales (Caissie 2006, Hauer et al. 2016) and are connected through a continuum of physicochemical processes across hierarchical spatial scales (Vannote et al. 1980). Studying the ecology of riverine species can be challenging, as the distribution of species in rivers can occur both independently of and connected to these dynamic processes (Fausch et al. 2002). Because of the complexity of river environments, multiscale sampling approaches are necessary to contextualize findings within small scales to the availability of habitats across larger scales (Fausch et al. 2002). Despite this, logistical constraints of sampling these ecosystems often means studies are conducted at intermediate scales and findings are extrapolated across unsampled units (Torgersen et al. 2022). Further, conservation management often occurs at mismatched scales to fisheries research, and few studies have effectively employed multiscale research in a way that can 150 be useful to conservation managers (Torgersen et al. 2022). These factors drove the conception of this dissertation. I studied the spatial and thermal ecology of Arctic grayling as a model species to gain insights about what drives (and over what scales) the spatiotemporal distribution of salmonids foraging in rivers. I consider these findings and their implications to improving conservation management in rivers. To achieve my objectives, I developed methods which have been useful in other systems but have been underutilized in river ecology. I then employed these methods at the reach, river, and watershed scales to demonstrate the importance of multiscale studies in rivers. Methodological developments In this dissertation I make several methodological contributions which can be useful to fish ecologists studying river systems. In my second chapter, I developed spatial capture-recapture (SCR) models which could be applied to large acoustic telemetry (AT) datasets in dendritic (branching) river networks. To-date, many of the methods for analysing AT data have been descriptive, focusing on where (and when) fish move (e.g. Niella et al. 2020; Flávio and Baktoft 2021). These approaches have been widely used in AT science but are limited in that they are unable to uncover the environmental drivers of the movements they characterize. Conversely, receiver-level regression and site occupancy modeling (e.g. O’Connor 2023) have been useful at quantifying these drivers but are limited to providing estimates of presence/absence and cannot account for the relative density of tagged animals in space. SCR modeling can incorporate both aspects, but to-date has been mostly limited to terrestrial systems and until this work had no clear guidance on how AT data could be integrated. SCR models use the full detection history of each tagged fish and can accommodate discontinuous receiver arrays, making them ideal for AT applications. To develop existing formulations for 151 application to my study system in the Parsnip watershed, two key developments had to be made. The first was the reliance of SCR models on Euclidean distance, which when applied within a dendritic watershed network would often predict activity centers on land. This was addressed by using Efford’s (2023) formulation of the model which replaces the distance term of the detection function with network distance. The second development was related to AT data specifically. Acoustic receivers do not have the ability to detect new, untagged individuals in their vicinity. This can bias results because estimates are dependent on the distribution of tagged fish, which may itself be biased based on the distribution of tagging effort. I addressed this by including a spatial layer of tagging effort, which allowed the detection function to thin predictions by the likelihood that an animal was already tagged at that location. Together, these developments presented a scalable framework which can be useful to ecologists in both river networks and 2D open-water systems. Further, this framework has flexible potential for use with remotely-sensed covariates. Over watershed scales, satellite sensing has proven useful for characterizing coarse habitat features (Glassic et al. 2024). Over river and reach scales (with potential at watershed scales), drones have emerged as a costeffective alternative to traditional aerial surveys (Dugdale 2016). For studies characterizing the physical and thermal habitats of freshwater fishes, recently developed integrated FLIR + RGB sensors offer the potential for collecting a variety of data in a single flight. However, the application of these combined data is still developing. In this dissertation, I make several contributions to drone work in rivers. For physical mapping, habitats can be delineated at reach scales by manually tracing orthomosaic images (Woodget et al. 2017). From a data processing standpoint, this is feasible at reach scales (2 km of habitats can be manually delineated in about a day). Extraction of exact habitat parameters (habitat areas, temperature summary statistics) using spatial algorithms 152 in R (e.g. Dunnington et al. 2023; Hijmans 2023) is computationally expensive, but doable within memory allowances at this scale. Extending the manual delineation approach to the river scale, while technically possible, presents analytical bottlenecks which can be prohibitive to projects that cover large areas. For example, scaling my drone surveys from the reach scale to the river scale required days instead of hours, but to manually delineate all of the habitats in the Anzac River, processing times would have scaled from hours/days to multiple weeks. The relative ease and affordability of drone methods for collecting large-scale imagery presents fisheries ecologists with abundant opportunities across large spatial scales, though methodologies for data processing over large scales are still lacking. I developed a pool index which could be used to survey orthomosaics to assign spatial point patterns to the distribution and characteristics of discrete pool habitats along the length of the river. In doing so, this reduced analysis times from multiple weeks to 2 days at the river scale, and the spatial points were readily exported and aggregated by Rkm for use in modeling. For studies in which discrete habitat features (e.g. pools, LWD, macrophyte patches) are of interest, this can be a useful approach for studying species across large scales in other systems. While this approach loses some resolution when compared to exact habitat delineation at the reach scale, it is feasible over large spatial extents and its usefulness was demonstrated in my third chapter. However, it is still a manual approach. Development of methods which can auto-delineate habitats from FLIR + RGB in a similar way to the NIR + RGB approach of Thompson et al. (2021) would be a useful step towards streamlining riverine drone work. Thermal habitat mapping in rivers is also still under development. FLIR applications from helicopters have been present for decades (Torgersen et al. 1999), but drone-payload FLIR cameras are newer. As such, their potential as a tool for fisheries ecologists is rapidly expanding. Most 153 applications to-date have focused on the identification of coldwater patches which represent discrete thermal refugia within the river (e.g. Casas-Mulet et al. 2020; Dugdale et al. 2015; Wilms and Whitworth 2016). However, applications of this method are still relatively few. The last comprehensive review of thermal drone techniques in rivers was published by Dugdale et al. (2019). The authors present both guidelines for data collection and cautions for interpretation of thermal data given the myriad factors that influence reflectivity in water (Dugdale et al. 2019). In my work, I used the FLIR rasters collected by the drone to represent the spatial distribution of temperatures at the reach scale (Chapter 3). Per existing recommendations, I used an instream temperature logger array to validate and correct my reach-scale thermal imagery, but a lack of a river-scale array precluded me from using the FLIR data in the river-scale models. Instream loggers have a high temporal resolution (5 minutes in my study) but limited spatial resolution (one point). FLIR imagery has a high spatial resolution (5 cm) but limited temporal resolution (one snapshot). In this dissertation, I experimented with a way to leverage the maximum resolution in each dimension to create reach-scale hourly rasters. Some early work in temporal interpolation of thermal imagery has been done in terrestrial systems (Weng et al. 2019; Peng et al. 2023), but to my knowledge has not been attempted in rivers. This assumed that the relative distribution of temperatures in flowing water remained similar throughout the day, and while it was not a perfect representation it presented a novel way to characterize the continuous data provided by the temperature loggers across the temporal snapshots provided by the FLIR rasters. Further research into the daily distribution of temperatures in lotic systems would aid in the further development of this method, and refinement could lead to potentially representing the daily thermal experience of riverine fishes. 154 Ecological insights and applications to conservation I used my combination of new and established methods to study a population of Arctic grayling of conservation concern in the Parsnip River watershed in north-central British Columbia, Canada. Following the impoundment of the Williston Reservoir in 1967, this fluvial population faced major declines due to extensive habitat loss and overfishing. Harvest regulations and long-term monitoring have led to some positive outcomes for the population; the number of Arctic grayling have doubled in the watershed over the past 25 years (Hagen and Stamford 2023). However, there are still data gaps that need to be filled before boots-on-the-ground conservation work can occur (e.g. coldwater engineering, land procurement for protection, riparian habitat restoration, spatial fishery closures; Stamford et al. 2017). Through my three data chapters, I identified eleven new insights about how Arctic grayling use their habitats in the Anzac River and the greater Parsnip River watershed. These insights and their implications for the conservation of this population are presented by chapter order in Table 5.1. The findings in this work can be grouped into two broad categories: habitat ecology and behavioural ecology. Table 5.1. Insights into Arctic grayling ecology found in this dissertation and their implications for the conservation of the Parsnip River watershed population. The umbrella term ‘Protection’ is used to describe conservation implications, but this can encompass habitat restoration, enhancement, land procurement for protection actions. Insight 1 2 Chapter Finding Conservation implication 2 Watershed-scale temperatures (11.1 - 17.1 °C) are important in the summer Protection or enhancement of coldwater-producing habitats across the watershed 2 Arctic grayling and bull trout occurrence are positively associated in the summer and winter, but not spring Protections of one species can benefit both 155 3 Pool habitats were important to Arctic grayling at both reach and river scales Protection of highlystructured instream habitats 3 Temperature (11.0 - 16.0 °C) was important at the reach scale, but uncertain at the river scale Protection of coldwaterproducing habitats; protection of highly-structured instream habitats 5 3 Upstream distance was important at the river scale, but uncertain at the reach scale Site fidelity of Arctic grayling may override river-scale temperature distributions, increased risk to river warming 6 4 Thermal preference of Arctic grayling was 10.1 - 13.0 °C Protection or enhancement of coldwater-producing habitats across the watershed 4 Heat transfer in Arctic grayling is slightly more efficient when cooling than warming Protection or enhancement of coldwater-producing habitats across the watershed 4 Heat transfer in Arctic grayling is more rapid in smaller fish Smaller, downstream fish, may expend greater energy reserves and show shorter excursion times beyond suitable thermal habitats 3 4 7 8 9 10 4 4 Arctic grayling use behavioural thermoregulation to maintain their body temperatures The effectiveness of behavioural thermoregulation was more strongly related to thermal heterogeneity in time (daily and seasonally) than in space Short-term resiliency to warming may come at tradeoff of long-term fitness of population Inseason fisheries management can exploit temperature-dependent fisheries closures over subdaily scales 156 11 4 Arctic grayling showed evidence of being onedirection (cooling) thermoregulators Protection or enhancement of coldwater-producing habitats Habitat ecology In this work, I investigated the thermal and physical habitat use of Arctic grayling during their summer feeding period which lasted from approximately July 1 – September 15 (with some variation between study years). I found positive associations with pool habitats at both the river and reach scales (Insight 3), though data on pool habitats were not available at the watershed scale. This finding corresponds well with prior work on Arctic grayling. Strong associations with pool habitats were described in Hughes (1992a, 1992b), McPhail (2007), Blackman (2002), and Zemlak and Langston (1998). While this finding was expected given what is known about Arctic grayling habitat associations, how these pool habitats interacted with temperature and upstream distance led to findings that were unexpected. Temperature was examined at all three scales of this study, and as the ‘master variable’ of fish ecology (Brett 1971) I expected to find a significant relationship across all scales. Temperature estimates from the thermal preference study (10.1 – 13.0 °C; Insight 7), the reach-scale study (11.0 – 16.0 °C; Insight 4), and the watershed-scale study (11.1 – 17.1 °C; Insight 1) all indicated that Arctic grayling were associated with cooler water temperatures and were in-line with estimates derived from other Arctic grayling studies (Hubert et al. 1985; Stewart et al. 2007; Hawkshaw 2011; Larocque et al. 2014; O’Connor 2023). However, I did not find a temperature signal at the river scale. Instead, I found that a combination of upstream distance and pool habitats were stronger predictors of Arctic grayling abundance (Insights 4 and 5). 157 Together, these findings had several implications to the conservation of this population. As evidenced at all scales except the river scale, the conservation of coldwater-producing habitats (e.g., headwater protection, coldwater engineering/augmentation, riparian planting) is important to Arctic grayling at both micro and macro scales in this watershed. Indeed, this approach is wellestablished and often favored by conservation practitioners working in river ecosystems (Kaylor et al. 2021; Kurylyk et al. 2015). In my third chapter, I found that upstream distance likely accounted for the effects of temperature as well as other unmeasured variables (e.g. forage density) that are also important to Arctic grayling. However, the river-scale findings draw attention to a further approach that can be employed to produce better outcomes for Arctic grayling than coldwater conservation alone. A similar finding to mine was identified in Chinook salmon in the John Day River basin in Oregon, where river-scale temperature distributions suggested that fish would have limited their habitat use to cold headwater reaches where stressful or lethal temperatures could be avoided (Torgersen et al. 1999). However, researchers found that thermal refugia associated with pool habitats, which allowed fish to escape high temperatures using discrete habitat features within otherwise unfavorable reaches, were driving habitat use at the river scale (Torgersen et al. 1999). Indeed, one of the core principles of multiscale riverscape sampling is that discrete habitat features (e.g. thermal refugia or groundwater upwelling) can have an overriding effect on the distribution of stream fishes (Baxter and Hauer 2000; Fausch et al. 2002). In my study, it is likely that the 1-Rkm resolution of the river-scale temperature models was too coarse to capture the overriding effects of thermal refugia associated with pools, while the pool metric itself was derived from the specific point patterns and characteristics of individual pools within the riverscape, leading to pools being a stronger predictor. 158 Coldwater conservation in rivers often focuses on the preservation of headwaters, but this finding highlights the importance of highly structured instream habitats as well (Kaylor et al. 2021). The other factor which was not identified in the Chinook salmon study (Torgersen et al. 1999), but was identified in my study, was upstream distance (Insight 5), and this may have specific implications to the conservation of Arctic grayling in the Anzac River. Fluvial Arctic grayling show high site fidelity, often returning to the same pool for feeding between years (Blackman 2002). Further, upstream distance is important to Arctic grayling dominance hierarchies (Hughes and Dill 1992; Hughes 1998), though its effects on abundance are uncertain (Fitzsimmons and Blackburn 2007; Fitzsimmons et al. 2009; Barker et al. 2011). While Arctic grayling will abandon their feeding sites if certain factors (e.g. turbidity) become limiting to their feeding success (Stamford et al. 2017), it is not certain to what degree or at what thresholds high temperatures may lead to this behavioural response. If site fidelity and maintaining their position within the river dominance hierarchy are more powerful drivers than seeking thermal refuge, Arctic grayling may be particularly susceptible to critical heat stress events. While this should be studied further with an explicit movement tracking study during a hot year (which may be possible using the 2021 acoustic telemetry data from Chapter 2), this provides a rationale to employ engineered coldwater augmentation (e.g., diversion of coldwater inputs, changes to channel structure) as a conservation strategy near highlyused feeding pools along the river (Kurylyk et al. 2015). It also provides rationale to employ management-based spatial fisheries closures. If Arctic grayling are highly reliant on discrete instream reaches, reducing angling pressures in high-traffic areas which overlap with highly structured instream habitats would provide better outcomes for the population over the long term. 159 Behavioural ecology While behavioural thresholds for abandoning feeding sites need further study, many of the insights from this dissertation provide further information on Arctic grayling conservation. I found that Arctic grayling used behavioural thermoregulation to regulate their body temperatures during their summer feeding period (Insight 9). My research provides novel evidence of Arctic grayling being single-direction thermoregulators (Insight 11), and that their bodies are marginally better at cooling than warming (Insight 7). As reported for other ectothermic species (e.g. Connoy et al. 2020), larger individuals took longer to warm after leaving favorable thermal habitats (Insight 8), which theory suggests would allow them to make longer excursions into warmer habitats and allocate their energy reserves more efficiently. Together, these findings provide new considerations for conservation prioritization (sensu Bayly et al. 2024; Hanson et al. 2019; Joseph et al. 2009). Arctic grayling are a slow-growing and slow-maturing species (McPhail 2007), and following upstream dominance hierarchy theory (Hughes and Dill 1992), younger, smaller fish in warm downstream reaches may be more susceptible to the effects of thermal stress (Insight 8). It is also important to consider that behavioural thermoregulation itself, while effective, may only offer short-term resilience to heat stress within generations (Kearney et al. 2009; Kelly 2019). Behavioural thermoregulation has been shown in other species to reduce population fitness across generations as avoiding exposure to high temperatures reduces the selective pressures of temperature on gene pools (Kearney et al. 2009; Kelly 2019). Conservation planning will have to consider these factors; smaller, younger fish may be less fit, but larger, more effective thermoregulators may supress population fitness. Depending on where population bottlenecks occur (which is not to my knowledge known for this population), conservation actions may be best applied from downstream to upstream (or vice-versa). It would be beneficial to evaluate this 160 population within the CEMPRA (Cumulative Effects Models for Prioritization of Restoration Actions; Bayly et al. 2024) framework. This would allow for various conservation actions (e.g. coldwater augmentation, riparian planting, enhancement) to be compared relative to their impacts at different life stages. In the specific Anzac River system, the findings related upstream distance support spatial management closures near the upper extent of the study area. As fish are distributed along an upstream gradient with the most fish upstream, this would reduce angling pressures on the largest and most fit individuals within the population. This approach would also protect areas of the river which are important staging areas in August for blue-listed bull trout ahead of their spawning migration. A spatial-seasonal closure of the areas around reach B would have the added benefit of reducing pressures on Arctic grayling during the hottest part of the summer, which circumvents some of the challenges with time-delays in existing temperature-dependent fishery closure strategies permitted under the British Columbia Drought and Water Scarcity Response Plan (BC government 2024). The combination of protections for both bull trout and Arctic grayling makes a compelling argument for this management action. On the topic of bull trout, the investigations in my study did not find strong evidence to support a predatory relationship. If I had, this would have been a complex issue to navigate as bull trout are protected at a higher order in British Columbia. Given my findings, protections of both species would have an additive effect. Beyond Arctic grayling While this discussion has so far largely focused on the ecological insights and conservation applications of my work to Arctic grayling, I consider aspects of this work to be important to the research of riverine species in general. I found that thermal heterogeneity in time was more important to thermoregulatory effectiveness than thermal heterogeneity in space, though this 161 finding may be an artifact of the scale and precision of my measurements of thermal variation within the spatial environment (Insight 10). While this was observed for Arctic grayling, it is likely that this pattern could be observed in other fluvial species. As mentioned above, it is possible this finding was influenced by the imperfect interpolated thermal raster data which I used to represent spatial heterogeneity. However, there are several factors which I consider germane to species thermoregulating in rivers. Much of thermoregulatory theory and its emphasis on spatial heterogeneity was developed using terrestrial ectotherms in arid environments where temperature differences between sun and cover patches can be quite pronounced over small spatial scales (e.g. Sears et al. 2016). Studies which have assessed thermoregulation in fishes have been less common, and those which have mostly investigated lakes and reservoir habitats with pronounced thermal stratification (e.g. Pepino et al. 2105, Armstrong et al. 2016; Biro 1998; Encina et al. 2008; Hitt et al. 2007; Amat-Trigo et al. 2023). Fewer have looked at rivers (Chiaramonte et al. 2016; Ritter et al. 2020), and these did not explicitly quantify the effectiveness of behavioural thermoregulation or spatial distribution of temperatures as in this study. While the spatial heterogeneity of thermal habitats in rivers is more structured than previously thought (Dzara et al. 2019), and groundwater upwelling can produce microhabitats which differ up to 10 °C from surrounding waters (Hare et al. 2021; Kanno et al. 2014; Power et al. 1999), river habitats are unique among freshwater environments in that they are also strongly heterogeneous along fine spatial and temporal axes (Caissie 2006). In my study reaches in the Anzac River, spatial heterogeneity of up to 4 °C was observed in the most pronounced cases, though 1 – 2 °C was typical. However, temporal heterogeneity could fluctuate nearly 8 °C over the diel cycle during the hottest part of the season, and discrete habitats could change up to 1 °C over just 90 minutes. 162 Responses and recovery from thermal stress in fishes is proportionate to both the magnitude and the duration of the stress (Farrell et al. 2008; Rezende et al. 2014). It follows then, that the proportion of the diel temperature cycle which falls below (and above) stress thresholds should be a key consideration when planning conservation action in rivers. If the proportion of recovery temperatures exceeds that of stress temperatures, then behavioural thermoregulation during the hot part of the day may be a sufficient strategy for mitigating thermal stress (provided critical temperature thresholds are not exceeded). Conversely, if the proportion of stress temperatures exceeds that of recovery temperatures, then intervention will be necessary. Management-based vectors to enact these interventions exist in the form of temperature-dependent fisheries closures (e.g. B.C. government 2015), though management timelines can operate over different scales than heat stress events. As such, temperature-based management strategies can be limited and the most effective solution to problems in thermal ecology will always be to slow down climate change. However, further research into species-specific stress thresholds and daily stress/recovery temperature ratios could set the groundwork for advocacy for policy-oriented conservation strategies (e.g. making temperature-based fisheries closures an existing part of angling regulations based on a watershed’s species). This approach, combined with the preservation of coldwater habitats, important instream habitats, a risk-based prioritization approach, and seasonal spatial fishery closures would form a robust conservation framework applicable to species in sensitive riverine systems. 163 164 References Abram, P. K., Boivin, G., Moiroux, J. & Brodeur. (2017). J. Behavioural effects of temperature on ectothermic animals: unifying thermal physiology and behavioural plasticity. Biol Rev 92, 1859– 1876. Alberta Government. 2018. Alberta Fish Sustainability Index (FSI) data provided by Fish & Wildlife Policy Branch. Alós, J., Aarestrup, K., Abecasis, D., Afonso, P., Alonso‐Fernandez, A., Aspillaga, E., Barcelo‐ Serra, M., Bolland, J., Cabanellas‐Reboredo, M., Lennox, R., McGill, R., Özgül, A., Reubens, J., & Villegas‐Ríos, D. (2022). Toward a decade of ocean science for sustainable development through acoustic animal tracking. Global Change Biology, 28(19), 5630–5653. https://doi.org/10.1111/gcb.16343 Amat‐Trigo, F., Andreou, D., Gillingham, P. K., & Britton, J. R. (2023). Behavioural thermoregulation in cold‐water freshwater fish: Innate resilience to climate warming? Fish and Fisheries, 24(1), 187–195. https://doi.org/10.1111/faf.12720 Anderson M. G., Walling D. E., and Bates P.D. (Editors) (1996). Floodplain Processes (Chichester, United Kingdom: Wiley). ISBN: 978-0-471-96679-1. Andrade, D. V., Gavira, R. S. B., & Tattersall, G. J. (2015). Thermogenesis in ectothermic vertebrates. Temperature, 2(4), 454–454. https://doi.org/10.1080/23328940.2015.1115570 Andrusak, G.F., Ward H. and R. Ptolemy. 2023. Provincial Fisheries Management: Drought Response Plan-2023. Prepared for the Ministry of Forests, Victoria, BC. Updated July 2023. 22 pp+ Angilletta, Jr., M. J., Niewiarowski, P. H., Dunham, A. E., Leaché, A. D., & Porter, W. P. (2004). Bergmann’s clines in ectotherms: Illustrating a life‐history perspective with sceloporine lizards. The American Naturalist, 164(6), E168–E183. https://doi.org/10.1086/425222 Angilletta, M. J. 2009. Thermal Adaptation: A Theoretical and Empirical Synthesis. Oxford University Press. Anttila, K., Casselman, M. T., Schulte, P. M., & Farrell, A. P. (2013). Optimum temperature in juvenile salmonids: Connecting subcellular indicators to tissue function and whole-organism thermal optimum. Physiological and Biochemical Zoology, 86(2), 245–256. https://doi.org/10.1086/669265 Armstrong, J. B., Fullerton, A. H., Jordan, C. E., Ebersole, J. L., Bellmore, J. R., Arismendi, I., Penaluna, B. E., & Reeves, G. H. (2021). The importance of warm habitat to the growth regime of cold-water fishes. Nature Climate Change, 11(4), 354–361. https://doi.org/10.1038/s41558-02100994-y Armstrong, J. B., Schindler, D. E., Ruff, C. P., Brooks, G. T., Bentley, K. E., & Torgersen, C. E. (2013). Diel horizontal migration in streams: Juvenile fish exploit spatial heterogeneity in thermal and trophic resources. Ecology, 94(9), 2066–2075. https://doi.org/10.1890/12-1200.1 B.C. Environmental Assessment Office (2023). Inspection record for Coastal GasLink Pipeline Project number 20230043_IR002. 165 Baccante, D. 2010. Further evidence of size gradients in Arctic grayling (Thymallus arcticus) along stream length. BC Journal of Ecosystems and Management 11(3):13–17. http://jem.forrex.org/index.php/jem/article/view/ Barker, O., N. Illar, M., and Foos, A. 2011. Peel Watershed Fish Habitat Assessment. Yukon Fish and Wildlife Branch Report TR-11-12, Whitehorse, Yukon, Canada Basson, C. H., Levy, O., Angilletta, M. J., & Clusella‐Trullas, S. (2017). Lizards paid a greater opportunity cost to thermoregulate in a less heterogeneous environment. Functional Ecology, 31(4), 856–865. https://doi.org/10.1111/1365-2435.12795 Baxter, C. V., & Hauer, F. R. (2000). Geomorphology, hyporheic exchange, and selection of spawning habitat by bull trout (Salvelinus confluentus). Canadian Journal of Fisheries and Aquatic Sciences, 57(7), 1470–1481. https://doi.org/10.1139/f00-056 Bayly, M., J., Tekatch, A., M., Rosenfeld, J., Jarvis, L., Enders, E., 2024. Cumulative Effects Model for Prioritizing Recovery Actions (CEMPRA). Documentation prepared by ESSA Technologies Ltd. for the BC Ministry of Environment. 75 + iv pp. Birtwell, I. K., Korstrom, J. S., & Brotherston, A. E. (2003). Laboratory studies on the effects of thermal change on the behaviour and distribution of juvenile chum salmon in sea water. Journal of Fish Biology, 62(1), 85–96. https://doi.org/10.1046/j.1095-8649.2003.00008.x Blackman, B. G. June 2001. A strategic plan for the conservation and restoration of Arctic grayling in the Williston Reservoir Watershed. Peace/Williston Fish and Wildlife Compensation Program, Report No. 241. 16pp plus appendices Blackman, B.G. 2002. Radio telemetry studies of Arctic grayling migrations to overwinter, spawning and summer feeding areas in the Parsnip River watershed. Peace/Williston Fish and Wildlife Compensation Program, Prince George, B.C. Bloomfield, E. J., Guzzo, M. M., Middel, T. A., Ridgway, M. S., & McMeans, B. C. (2022). Seasonality can affect ecological interactions between fishes of different thermal guilds. Frontiers in Ecology and Evolution, 10, 986459. https://doi.org/10.3389/fevo.2022.986459 Blouin-Demers, G., & Nadeau, P. (2005). The cost-benefit model of thermoregulation does not predict lizard thermoregulatory behavior. Ecology, 86(3), 560–566. https://www.jstor.org/stable/3450648 Blouin-Demers, G., & Weatherhead, P. J. (2001). Thermal Ecology of Black Rat Snakes (Elaphe obsoleta) in a Thermally Challenging Environment. Ecology, 82(11), 3025–3043. https://doi.org/10.1890/0012-9658(2001)082[3025:TEOBRS]2.0.CO;2. Blouin-Demers, G., Bjorgan, L. P. G., & Weatherhead, P. J. (2007). Changes in habitat use and movement patterns with body size in black ratsnakes (Elaphe obsoleta). Herpetologica, 63(4), 421–429. https://doi.org/10.1655/0018-0831(2007)63[421:CIHUAM]2.0.CO;2 Boehrer, B., and M. Schultze (2008), Stratification of lakes, Rev. Geophys., 46, RG2005, doi:10.1029/2006RG000210. 166 Borchers, D. L., & Efford, M. G. (2008). Spatially explicit maximum likelihood methods for capture-recapture studies. Biometrics, 64(2), 377–385. https://doi.org/10.1111/j.15410420.2007.00927.x Brennan, S. R., Schindler, D. E., Cline, T. J., Walsworth, T. E., Buck, G., & Fernandez, D. P. (2019). Shifting habitat mosaics and fish production across river basins. Science, 364(6442), 783–786. https://doi.org/10.1126/science.aav4313 Brewitt, K. S., Danner, E. M., & Moore, J. W. (2017). Hot eats and cool creeks: Juvenile Pacific salmonids use mainstem prey while in thermal refuges. Canadian Journal of Fisheries and Aquatic Sciences, 74, 1588–1602. https://doi.org/10.1139/cjfas-2016-0395 Brewster, C. L., Sikes, R. S., & Gifford, M. E. (2013). Quantifying the cost of thermoregulation: Thermal and energetic constraints on growth rates in hatchling lizards. Functional Ecology, 27(2), 490–497. https://doi.org/10.1111/1365-2435.12066 Brooks, M. E., Kasper Kristensen, Koen J. van Benthem, Arni Magnusson, Casper W. Berg, Anders Nielsen, Hans J. Skaug, Martin Maechler and Benjamin M. Bolker (2017). glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling. The R Journal, 9(2), 378-400. doi:10.32614/RJ-2017-066. Brosset, P., Averty, A., Mathieu-Resuge, M., Schull, Q., Soudant, P., & Lebigre, C. (2023). Fish morphometric body condition indices reflect energy reserves but other physiological processes matter. Ecological Indicators, 154, 110860. https://doi.org/10.1016/j.ecolind.2023.110860 Brown, G. P., & Weatherhead, P. J. (2000). Thermal ecology and sexual size dimorphism in northern water snakes, Nerodia sipedon. Ecological Monographs, 70(2), 311–330. https://doi.org/10.2307/2657179 Brownscombe, J. W., Griffin, L. P., Brooks, J. L., Danylchuk, A. J., Cooke, S. J., & Midwood, J. D. (2022). Applications of telemetry to fish habitat science and management. Canadian Journal of Fisheries and Aquatic Sciences, 79(8), 1347–1359. https://doi.org/10.1139/cjfas-2021-0101 Brett, J. R. (1971). Energetic responses of salmon to temperature. A study of some thermal relations in the physiology and freshwater ecology of sockeye salmon (Oncorhynchus nerka). American Zoologist, 11(1), 99–113. https://doi.org/10.1093/icb/11.1.99 Bunnell, F., D.F. Fraser and A. Harcombe. 2009. Increasing effectiveness of conservation decisions with existing data: a simple system and its application. Natural Areas Journal. Volume 29 (1) 7990. Burnham, K. P., Anderson, D. R., & Burnham, K. P. (2002). Model selection and multimodel inference: A practical information-theoretic approach (2nd ed). Springer. Busch, S., Kirillin, G., & Mehner, T. (2012). Plasticity in habitat use determines metabolic response of fish to global warming in stratified lakes. Oecologia, 170(1), 275–287. https://doi.org/10.1007/s00442-012-2286-z Buzby, K. M., & Deegan, L. A. (2000). Inter-annual fidelity to summer feeding sites in arctic grayling. Environmental Biology of Fishes, 59(3), 319–327. https://doi.org/10.1023/A:1007626507936 167 Cahill (2015). Alberta Environment and Parks and Alberta Conservation Association. 2015. Status of the Arctic Grayling (Thymallus arcticus) in Alberta: Update 2015. Alberta Environment and Parks. Alberta Wildlife Status Report No. 57 (Update 2015). Edmonton, AB. 96 pp. Caissie, D. (2006). The thermal regime of rivers: A review. Freshwater Biology, 51(8), 1389–1406. https://doi.org/10.1111/j.1365-2427.2006.01597.x Carrillo‐Longoria, J., Gaylord, G., Andrews, L., & Powell, M. (2024). Effect of temperature on growth, survival, and chronic stress responses of Arctic Grayling juveniles. Transactions of the American Fisheries Society, 153(1), 3–22. https://doi.org/10.1002/tafs.10453 Casas-Mulet, R., Pander, J., Ryu, D., Stewardson, M. J., & Geist, J. (2020). Unmanned aerial vehicle (Uav)-based thermal infra-red (Tir) and optical imagery reveals multi-spatial scale controls of cold-water areas over a groundwater-dominated riverscape. Frontiers in Environmental Science, 8, 64. https://doi.org/10.3389/fenvs.2020.00064 Chapman, B. B., Hulthén, K., Brodersen, J., Nilsson, P. A., Skov, C., Hansson, L.-A., & Brönmark, C. (2012). Partial migration in fishes: causes and consequences. Journal of Fish Biology, 81(2), 456–478. https://doi.org/10.1111/j.1095-8649.2012.03342.x Chrétien, L., Théau, J., & Ménard, P. (2016). Visible and thermal infrared remote sensing for the detection of white‐tailed deer using an unmanned aerial system. Wildlife Society Bulletin, 40(1), 181–191. https://doi.org/10.1002/wsb.629 Christensen, E. A. F., Andersen, L. E. J., Bergsson, H., Steffensen, J. F., & Killen, S. S. (2021). Shuttle-box systems for studying preferred environmental ranges by aquatic animals. Conservation Physiology, 9(1), coab028. https://doi.org/10.1093/conphys/coab028 Christian, K. A., & Weavers, B. W. (1996). Thermoregulation of monitor lizards in Australia: An evaluation of methods in thermal biology. Ecological Monographs, 66(2), 139–157. https://doi.org/10.2307/2963472 Clarke, A. D., Telmer, K. H., & Shrimpton, J. M. (2007). Habitat use and movement patterns for a fluvial species, the Arctic grayling, in a watershed impacted by a large reservoir: Evidence from otolith microchemistry. Journal of Applied Ecology, 44(6), 1156–1165. https://doi.org/10.1111/j.1365-2664.2007.01350.x Clements, W. H., Herbst, D. B., Hornberger, M. I., Mebane, C. A., & Short, T. M. (2021). Longterm monitoring reveals convergent patterns of recovery from mining contamination across 4 western US watersheds. Freshwater Science, 40(2), 407–426. https://doi.org/10.1086/714575 Comte, L., & Grenouillet, G. (2013). Do stream fish track climate change? Assessing distribution shifts in recent decades. Ecography, 36(11), 1236–1246. https://doi.org/10.1111/j.16000587.2013.00282.x Connoy, J. W. H., Leivesley, J. A., Brooks, R. J., Litzgus, J. D., & Rollinson, N. (2020). Body size of ectotherms constrains thermal requirements for reproductive activity in seasonal environments. Canadian Journal of Zoology, 98(10), 651–659. https://doi.org/10.1139/cjz-2019-0254 Crook, D. A. (2004). A method for externally attaching radio transmitters to minimize dermal irritation. Journal of Fish Biology, 64(1), 258–261. https://doi.org/10.1111/j.10958649.2004.00282.x 168 Crossin, G. T., Heupel, M. R., Holbrook, C. M., Hussey, N. E., Lowerre‐Barbieri, S. K., Nguyen, V. M., Raby, G. D., & Cooke, S. J. (2017). Acoustic telemetry and fisheries management. Ecological Applications, 27(4), 1031–1049. https://doi.org/10.1002/eap.1533 Dahlke, F. T., Wohlrab, S., Butzin, M., & Pörtner, H.-O. (2020). Thermal bottlenecks in the life cycle define climate vulnerability of fish. Science, 369(6499), 65–70. https://doi.org/10.1126/science.aaz3658 Dahms, C., & Killen, S. (2022). Temperature change effects on marine fish range shifts: A metaanalysis of ecological and methodological predictors [Preprint]. Life Sciences. https://doi.org/10.32942/X2C88T Díaz, F., Re, A. D., González, R. A., Sánchez, L. N., Leyva, G., & Valenzuela, F. (2007). Temperature preference and oxygen consumption of the largemouth bass Micropterus salmoides (Lacépède) acclimated to different temperatures: Temperature preference and oxygen consumption of the largemouth bass. Aquaculture Research, 38(13), 1387–1394. https://doi.org/10.1111/j.13652109.2007.01817.x Donaldson, M. R., Cooke, S. J., Patterson, D. A., & Macdonald, J. S. (2008). Cold shock and fish. Journal of Fish Biology, 73(7), 1491–1530. https://doi.org/10.1111/j.1095-8649.2008.02061.x Donaldson, M. R., Cooke, S. J., Patterson, D. A., Hinch, S. G., Robichaud, D., Hanson, K. C., Olsson, I., Crossin, G. T., English, K. K., & Farrell, A. P. (2009). Limited behavioural thermoregulation by adult upriver-migrating sockeye salmon (Oncorhynchus nerka) in the Lower Fraser River, British Columbia. Canadian Journal of Zoology, 87(6), 480–490. https://doi.org/10.1139/Z09-032 Dorazio, R. M., & Karanth, K. U. (2017). A hierarchical model for estimating the spatial distribution and abundance of animals detected by continuous-time recorders. PLOS ONE, 12(5), e0176966. https://doi.org/10.1371/journal.pone.0176966 Dressler, T.L., Han Lee, V., Klose, K. et al. Thermal tolerance and vulnerability to warming differ between populations of wild Oncorhynchus mykiss near the species’ southern range limit. Sci Rep 13, 14538 (2023). https://doi.org/10.1038/s41598-023-41173-7 Drusch, M., Wood, E. F., and Gao, H. (2005). Observation operators for the direct assimilation of TRMM microwave imager retrieved soil moisture. Geophys. Res. Lett. 32, 32–35. doi: 10.1029/2005GL023623 Dudgeon, D., Arthington, A. H., Gessner, M. O., Kawabata, Z., Knowler, D. J., Lévêque, C., Naiman, R. J., Prieur‐Richard, A., Soto, D., Stiassny, M. L. J., & Sullivan, C. A. (2006). Freshwater biodiversity: Importance, threats, status and conservation challenges. Biological Reviews, 81(2), 163–182. https://doi.org/10.1017/S1464793105006950 Dugdale, S. J., Bergeron, N. E., & St-Hilaire, A. (2015). Spatial distribution of thermal refuges analysed in relation to riverscape hydromorphology using airborne thermal infrared imagery. Remote Sensing of Environment, 160, 43–55. https://doiorg.prxy.lib.unbc.ca/10.1016/j.rse.2014.12.021 169 Dugdale, S. J. (2016). A practitioner's guide to thermal infrared remote sensing of rivers and streams: Recent advances, precautions and considerations. Wiley Interdisciplinary Reviews: Water, 3, 251–268. https://doi-org.prxy.lib.unbc.ca/10.1002/wat2.1135 Dugdale, S. J., Kelleher, C. A., Malcolm, I. A., Caldwell, S., & Hannah, D. M. (2019). Assessing the potential of drone‐based thermal infrared imagery for quantifying river temperature heterogeneity. Hydrological Processes, 33(7), 1152–1163. https://doi.org/10.1002/hyp.13395 Dzara, J. R., Neilson, B. T., & Null, S. E. (2019). Quantifying thermal refugia connectivity by combining temperature modeling, distributed temperature sensing, and thermal infrared imaging. Hydrology and Earth System Sciences, 23(7), 2965–2982. https://doi.org/10.5194/hess-23-29652019 Efford, M. G., Boulanger, J., Spencer, N. J., Warburton, B. and Stenhouse, G. In prep. Joint analysis of telemetry and capture–recapture data. Efford, M. G. (2023). secrlinear: Spatially Explicit Capture--Recapture for Linear Habitats. R package version 1.1.4. https://CRAN.R-project.org/package=secrlinear/ Efford, M. G. (2023). secr: Spatially explicit capture-recapture models. R package version 4.5.8. https://CRAN.R-project.org/package=secr Efford, M. G., Borchers, D. L., & Byrom, A. E. (2009). Density estimation by spatially explicit capture–recapture: Likelihood-based methods. In D. L. Thomson, E. G. Cooch, & M. J. Conroy (Eds.), Modeling Demographic Processes In Marked Populations (pp. 255–269). Springer US. https://doi.org/10.1007/978-0-387-78151-8_11 Eliason, E. J., Clark, T. D., Hague, M. J., Hanson, L. M., Gallagher, Z. S., Jeffries, K. M., Gale, M. K., Patterson, D. A., Hinch, S. G., & Farrell, A. P. (2011). Differences in thermal tolerance among sockeye salmon populations. Science, 332(6025), 109–112. https://doi.org/10.1126/science.1199158 Elsner, R. A., & Shrimpton, J. M. (2019). Behavioural changes during the parr–smolt transformation in coho salmon Oncorhynchus kisutch: is it better to be cool? Journal of Fish Biology, jfb.14069. https://doi.org/10.1111/jfb.14069 Enrique Ángeles-González, L., Martínez-Meyer, E., Yañez-Arenas, C., Velázquez-Abunader, I., Garcia-Rueda, A., Díaz, F., Tremblay, N., Antonio Flores-Rivero, M., Gebauer, P., & Rosas, C. (2020). Using realized thermal niche to validate thermal preferences from laboratory studies. How do they stand? Ecological Indicators, 118, 106741. https://doi.org/10.1016/j.ecolind.2020.106741 Everett, K., C. Cadrin and T. Lea. 2009. Ecosystem prioritization technical methods paper. Ministry of Environment Technical Paper. Farrell, A. P. et al. Pacific Salmon in Hot Water: Applying Aerobic Scope Models and Biotelemetry to Predict the Success of Spawning Migrations. Physiol Biochem Zool 81, 697–709 (2008). Fausch, K. D., Torgersen, C. E., Baxter, C. V., & Li, H. W. (2002). Landscapes to Riverscapes: Bridging the Gap between Research and Conservation of Stream Fishes. BioScience, 52(6), 483. https://doi.org/10.1641/0006-3568(2002)052[0483:LTRBTG]2.0.CO;2 170 Fechhelm, R. G., & Neill, W. H. (1982). Predicting body-core temperature in fish subjected to fluctuating ambient temperature. Physiological Zoology, 55(3), 229–239. https://doi.org/10.1086/physzool.55.3.30157887 Fieberg, J., Signer, J., Smith, B., & Avgar, T. (2021). A ‘How to’ guide for interpreting parameters in habitat‐selection analyses. Journal of Animal Ecology, 90(5), 1027–1043. https://doi.org/10.1111/1365-2656.13441 Fisheries and Oceans Canada. 2020. Recovery Strategy for the Bull Trout (Salvelinus confluentus), Saskatchewan-Nelson Rivers populations, in Canada [Final]. Species at Risk Act Recovery Strategy Series. Fisheries and Oceans Canada, Ottawa. viii + 130 pp. Flávio, H. and Baktoft, H. (2021) actel: Standardised analysis of AT data from animals moving through receiver arrays. Methods in Ecology and evolution. doi: 10.1111/2041-210X.13503 Frechette, D. M., Dugdale, S. J., Dodson, J. J., & Bergeron, N. E. (2018). Understanding summertime thermal refuge use by adult Atlantic salmon using remote sensing, river temperature monitoring, and acoustic telemetry. Canadian Journal of Fisheries and Aquatic Sciences, 75, 1999– 2010. https://doi.org/10.1139/cjfas-2017–0422 Fullerton, A. H., Torgersen, C. E., Lawler, J. J., Faux, R. N., Steel, E. A., Beechie, T. J., Ebersole, J. L., & Leibowitz, S. G. (2015). Rethinking the longitudinal stream temperature paradigm: Region-wide comparison of thermal infrared imagery reveals unexpected complexity of river temperatures: COMPLEXITY IN LONGITUDINAL THERMAL PROFILES WITHIN AND AMONG RIVERS. Hydrological Processes, 29(22), 4719–4737. https://doi.org/10.1002/hyp.10506 Fullerton, A. H., Torgersen, C. E., Lawler, J. J., Steel, E. A., Ebersole, J. L., & Lee, S. Y. (2018). Longitudinal thermal heterogeneity in rivers and refugia for coldwater species: Effects of scale and climate change. Aquatic Sciences, 80(1), 3. https://doi.org/10.1007/s00027-017-0557-9 Gelman, A., & Hill, J. (2006). Data analysis using regression and multilevel/hierarchical models (1st ed.). Cambridge University Press. https://doi.org/10.1017/CBO9780511790942 Glassic, H. C., McGwire, K. C., Macfarlane, W. W., Rasmussen, C., Bouwes, N., Wheaton, J. M., & Al‐Chokhachy, R. (2024). From pixels to riverscapes: How remote sensing and geospatial tools can prioritize riverscape restoration at multiple scales. WIREs Water, 11(3), e1716. https://doi.org/10.1002/wat2.1716 Goniea, T. M., Keefer, M. L., Bjornn, T. C., Peery, C. A., Bennett, D. H., & Stuehrenberg, L. C. (2006). Behavioral thermoregulation and slowed migration by adult fall Chinook salmon in response to high Columbia River water temperatures. Transactions of the American Fisheries Society, 135, 408–419. https://doi.org/10.1577/T04-113.1 Goyer, K., Bertolo, A., Pépino, M., & Magnan, P. (2014). Effects of lake warming on behavioural thermoregulatory tactics in a cold-water stenothermic fish. PLoS ONE, 9(3), e92514. https://doi.org/10.1371/journal.pone.0092514 Guerrero, A. M., McAllister, R. R. J., Corcoran, J., & Wilson, K. A. (2013). Scale mismatches, conservation planning, and the value of social‐network analyses. Conservation Biology, 27(1), 35– 44. https://doi.org/10.1111/j.1523-1739.2012.01964.x 171 Gunderson, A. R., and M. Leal. 2016. A conceptual framework for understanding thermal constraints on ectotherm activity with implications for predicting responses to global change. Ecology Letters 19:111 120. Guzzo, M. M., & Blanchfield, P. J. (2017). Climate change alters the quantity and phenology of habitat for lake trout (Salvelinus namaycush) in small Boreal Shield lakes. Canadian Journal of Fisheries and Aquatic Sciences, 74(6), 871–884. https://doi.org/10.1139/cjfas-2016-0190 Haesemeyer, M. (2020). Thermoregulation in fish. Molecular and Cellular Endocrinology, 518, 110986. https://doi.org/10.1016/j.mce.2020.110986 Hagen, J. & Gantner, N. (2019) Abundance and trend of Arctic grayling (Thymallus arcticus) in index sites of the Parsnip River watershed, 1995-2019. FWCP Report. Project No. PEA-F20-2959 Hagen, J., and Stamford, M. (2017). FWCP Arctic grayling monitoring framework for the Williston Reservoir Watershed. Report prepared for the Fish and Wildlife Compensation Program – Peace Region Hagen, J., Spendlow, I., & Pillipow, R. (2020). Critical spawning habitats and abundance of Bull trout in the Williston Reservoir watershed, 2019 (PEA-F20-F-2956). Hagen, J., and S. Weber. 2019. Limiting Factors, Enhancement Potential, Critical Habitats, and Conservation Status for Bull Trout of the Williston Reservoir Watershed: Information Synthesis and Recommended Monitoring Framework. Report prepared for the Fish and Wildlife Compensation Program – Peace Region, Prince George, BC. Hagen, J., and M. Stamford. 2023. Parsnip Arctic Grayling Critical Habitats and Abundance: 20182022 Final Report. Report prepared for the Fish and Wildlife Compensation Program –Peace Region, Prince George, BC. FWCP Project No. PEA-F23-F-3631 Hartig F (2022). _DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models_. R package version 0.4.6, . Hauer, F. R., Locke, H., Dreitz, V. J., Hebblewhite, M., Lowe, W. H., Muhlfeld, C. C., Nelson, C. R., Proctor, M. F., & Rood, S. B. (2016). Gravel-bed river floodplains are the ecological nexus of glaciated mountain landscapes. Science Advances, 2(6), e1600026. https://doi.org/10.1126/sciadv.1600026 Hawkshaw, S.C.F. and J.M. Shrimpton. 2014. Temperature preference and distribution of juvenile Arctic grayling (Thymallus arcticus) in the Williston Watershed, British Columbia, Canada. Fish and Wildlife Compensation Program – Peace Region Report No. 366. 40 pp plus appendices Haydt, N. T., Hocking, D. J., & Sterrett, S. (2022). Spatial capture–recapture derived turtle capture probabilities and densities in the Chesapeake and Ohio canal. Journal of Herpetology, 56(2). https://doi.org/10.1670/21-026 Hesselbarth, M.H.K., Sciaini, M., With, K.A., Wiegand, K., Nowosad, J. 2019. landscapemetrics: an open-source R tool to calculate landscape metrics. - Ecography 42:1648-1657 (v0.0). Hijmans R (2023). _terra: Spatial Data Analysis_. R package version 1.7-3, . 172 Hitt, N. P., Snook, E. L., & Massie, D. L. (2017). Brook trout use of thermal refugia and foraging habitat influenced by brown trout. Canadian Journal of Fisheries and Aquatic Sciences, 74, 406– 418. https://doi.org/10.1139/cjfas-2016-0255 Hoenner, X., Huveneers, C., Steckenreuter, A., Simpfendorfer, C., Tattersall, K., Jaine, F., Atkins, N., Babcock, R., Brodie, S., Burgess, J., Campbell, H., Heupel, M., Pasquer, B., Proctor, R., Taylor, M. D., Udyawer, V., & Harcourt, R. (2018). Australia’s continental-scale acoustic tracking database and its automated quality control process. Scientific Data, 5(1), 170206. https://doi.org/10.1038/sdata.2017.206 Hodge, B. W., Battige, K. D., & Rogers, K. B. (2017). Seasonal and temperature‐related movement of Colorado River cutthroat trout in a low‐elevation, Rocky Mountain stream. Ecology and Evolution, 7(7), 2346–2356. https://doi.org/10.1002/ece3.2847 Hoover, S. E. R., Shannon, L. G. W., & Ackerman, J. D. (2007). The effect of riparian condition on invertebrate drift in mountain streams. Aquatic Sciences, 69(4), 544–553. https://doi.org/10.1007/s00027-007-0909-y Hostetter, N. J., & Royle, J. A. (2020). Movement-assisted localization from acoustic telemetry data. Movement Ecology, 8(1), 15. https://doi.org/10.1186/s40462-020-00199-6 Howell, P. J., Dunham, J. B., & Sankovich, P. M. (2010). Relationships between water temperatures and upstream migration, cold water refuge use, and spawning of adult bull trout from the Lostine River, Oregon, USA. Ecology of Freshwater Fish, 19(1), 96–106. https://doi.org/10.1111/j.1600-0633.2009.00393.x Hubert, W. A., Helzner, R. S., Lee, L. A., & Nelson, P. C. (1985). Habitat suitability index models and instream flow suitability curves: Arctic grayling riverine populations (82/10.110). U.S. Fish and Wildlife Service. https://pubs.usgs.gov/publication/fwsobs82_10_110 Huey, R. B., & Slatkin, M. (1976). Cost and benefits of lizard thermoregulation. The Quarterly Review of Biology, 51(3), 363–384. https://doi.org/10.1086/409470 Huey, R. B., Hertz, P. E., & Sinervo, B. (2003). Behavioral drive versus behavioral inertia in evolution: A null model approach. The American Naturalist, 161(3), 357–366. https://doi.org/10.1086/346135 Huff, D.D., Hubler, S.L., and Borisenko, A.N. 2005. Using Field Data to Estimate the Realized Thermal Niche of Aquatic Vertebrates. North American Journal of Fisheries Management 25(1): 346–360. doi:10.1577/m03-231.1 Huey, R. B., Kearney, M. R., Krockenberger, A., Holtum, J. A. M., Jess, M., & Williams, S. E. (2012). Predicting organismal vulnerability to climate warming: Roles of behaviour, physiology and adaptation. Philosophical Transactions of the Royal Society B: Biological Sciences, 367(1596), 1665–1679. https://doi.org/10.1098/rstb.2012.0005 Hughes, N. F. (1992a). Ranking of Feeding Positions by Drift-Feeding Arctic Grayling (Thymallus arcticus) in Dominance Hierarchies. Canadian Journal of Fisheries and Aquatic Sciences, 49(10), 1994–1998. https://doi.org/10.1139/f92-222 Hughes, N. F. (1992b). Selection of Positions by Drift-Feeding Salmonids in Dominance Hierarchies: Model and Test for Arctic Grayling (Thymallus arcticus) in Subarctic Mountain 173 Streams, Interior Alaska. Canadian Journal of Fisheries and Aquatic Sciences, 49(10), 1999–2008. https://doi.org/10.1139/f92-223 Hughes, N. F., & Reynolds, J. B. (1994). Why do Arctic grayling (Thymallus arcticus) get bigger as you go upstream? Canadian Journal of Fisheries and Aquatic Sciences, 51(10), 2154–2163. https://doi.org/10.1139/f94-216 Hughes, N. F. (1998). A Model of Habitat Selection by Drift-Feeding Stream Salmonids at Different Scales. Ecology, 79(1), 281–294. https://doi.org/10.1890/00129658(1998)079[0281:AMOHSB]2.0.CO;2 Hughes, N. F. (1999). Population processes responsible for larger-fish-upstream distribution patterns of Arctic grayling (Thymallus arcticus) in interior Alaskan runoff rivers. Canadian Journal of Fisheries and Aquatic Sciences, 56(12), 2292–2299. https://doi.org/10.1139/f99-157 Hughes, N. F. (2000). Testing the ability of habitat selection theory to predict interannual movement patterns of a drift‐feeding salmonid. Ecology of Freshwater Fish, 9(1–2), 4–8. https://doi.org/10.1034/j.1600-0633.2000.90102.x Hughes, N. F., & Dill, L. M. (1990). Position choice by drift-feeding salmonids: Model and test for arctic grayling (Thymallus arcticus) in subarctic mountain streams, interior Alaska. Canadian Journal of Fisheries and Aquatic Sciences, 47(10), 2039–2048. https://doi.org/10.1139/f90-228 Hussey, N. E., Kessel, S. T., Aarestrup, K., Cooke, S. J., Cowley, P. D., Fisk, A. T., Harcourt, R. G., Holland, K. N., Iverson, S. J., Kocik, J. F., Mills Flemming, J. E., & Whoriskey, F. G. (2015). Aquatic animal telemetry: A panoramic window into the underwater world. Science, 348(6240), 1255642. https://doi.org/10.1126/science.1255642 Isaak, D. J., & Young, M. K. (2023). Cold-water habitats, climate refugia, and their utility for conserving salmonid fishes. Canadian Journal of Fisheries and Aquatic Sciences, 80(7), 1187– 1206. https://doi.org/10.1139/cjfas-2022-0302 Iverson, J. B., Higgins, H., Sirulnik, A., & Griffiths, C. (1997). Local and geographic variation in the reproductive biology of the snapping turtle (Chelydra serpentina). Herpetologica, 53(1), 96– 117. https://www.jstor.org/stable/3893247 Johnson, Douglas H., "The Comparison of Usage and Availability Measurements for Evaluating Resource Preference" (1980). USGS Northern Prairie Wildlife Research Center. 198. https://digitalcommons.unl.edu/usgsnpwrc/198 Jonsson, B. (2023). Thermal effects on ecological traits of salmonids. Fishes, 8(7), 337. https://doi.org/10.3390/fishes8070337 Kanno, Y., Harris, A. C., Kishida, O., Utsumi, S., & Uno, H. (2022). Complex effects of body length and condition on within‐tributary movement and emigration in stream salmonids. Ecology of Freshwater Fish, 31(2), 317–329. https://doi.org/10.1111/eff.12632 Karameta, E., Gavriilidi, I., Sfenthourakis, S., & Pafilis, P. (2023). Seasonal variation in the thermoregulation pattern of an insular agamid lizard. Animals, 13(20), 3195. https://doi.org/10.3390/ani13203195 174 Kays, R., Sheppard, J., Mclean, K., Welch, C., Paunescu, C., Wang, V., Kravit, G., & Crofoot, M. (2019). Hot monkey, cold reality: Surveying rainforest canopy mammals using drone-mounted thermal infrared sensors. International Journal of Remote Sensing, 40(2), 407–419. https://doi.org/10.1080/01431161.2018.1523580 Kearney, M., Shine, R., & Porter, W. P. (2009). The potential for behavioral thermoregulation to buffer “cold-blooded” animals against climate warming. Proceedings of the National Academy of Sciences, 106(10), 3835–3840. https://doi.org/10.1073/pnas.0808913106 Keefer, M. L., Clabough, T. S., Jepson, M. A., Johnson, E. L., Peery, C. A., & Caudill, C. C. (2018). Thermal exposure of adult Chinook salmon and steelhead: Diverse behavioral strategies in a large and warming river system. PLoS One, 13, e0204274. https://doi.org/10.1371/journal.pone.020474 Kelly, M. (2019). Adaptation to climate change through genetic accommodation and assimilation of plastic phenotypes. Philosophical Transactions of the Royal Society B: Biological Sciences, 374(1768), 20180176. https://doi.org/10.1098/rstb.2018.0176 Kelly, B. B., Siepker, M. J., & Weber, M. J. (2023). Effects of non-native Salmo trutta and multiscale habitat factors on native fishes in the Driftless Area. Canadian Journal of Fisheries and Aquatic Sciences, cjfas-2022-0213. https://doi.org/10.1139/cjfas-2022-0213 Kennedy, T. A., Yackulic, C. B., Cross, W. F., Grams, P. E., Yard, M. D., & Copp, A. J. (2014). The relation between invertebrate drift and two primary controls, discharge and benthic densities, in a large regulated river. Freshwater Biology, 59(3), 557–572. https://doi.org/10.1111/fwb.12285 Kurylyk, B. L., MacQuarrie, K. T. B., Linnansaari, T., Cunjak, R. A., & Curry, R. A. (2015). Preserving, augmenting, and creating cold‐water thermal refugia in rivers: Concepts derived from research on the Miramichi River, New Brunswick (Canada). Ecohydrology, 8(6), 1095–1108. https://doi.org/10.1002/eco.1566 Langston, A.R., and Cubberley, J.C. 2008. Assessing the origin of bull trout spawners in the Misinchinka River and the river’s potential as a redd count index system: 2004 and 2005 radiotelemetry results. Peace/Williston Fish and Wildlife Compensation Program Report No. 317 Larocque, S.M., Hatry, C., & Enders, E.C. 2014. Development of habitat suitability indices and bioenergetics models for Arctic grayling (Thymallus arcticus). Canadian technical report of fisheries and aquatic sciences 3097. Larsen, H. L., Møller-Lassesen, K., Enevoldsen, E. M. E., Madsen, S. B., Obsen, M. T., Povlsen, P., Bruhn, D., Pertoldi, C., & Pagh, S. (2023). Drone with mounted thermal infrared cameras for monitoring terrestrial mammals. Drones, 7(11), 680. https://doi.org/10.3390/drones7110680 Larsson, S. Thermal preference of Arctic charr, Salvelinus alpinus, and brown trout, Salmo trutta – implications for their niche segregation. Environ Biol Fish 73, 89–96 (2005). https://doi.org/10.1007/s10641-004-5353-4 Lennox, R. J., Chapman, J. M., Twardek, W. M., Broell, F., Bøe, K., Whoriskey, F. G., Fleming, I. A., Robertson, M., & Cooke, S. J. (2019). Biologging in combination with biotelemetry reveals behavior of Atlantic salmon following exposure to capture and handling stressors. Canadian Journal of Fisheries and Aquatic Sciences, 76(12), 2176–2183. https://doi.org/10.1139/cjfas2018-0477 175 Lennox, R.J., Eldøy, S.H., Dahlmo, L.S. et al. Acoustic accelerometer transmitters and their growing relevance to aquatic science. Mov Ecol 11, 45 (2023). https://doi.org/10.1186/s40462023-00403-3 Levin S.A. (1992). The problem of pattern and scale in ecology. Ecology 73: 1943–1967. Liuzzo, M., Borella, S., Ottonello, D., Arizza, V., & Malavasi, S. (2021). Population abundance, structure and movements of the European pond turtle, Emys orbicularis (Linnaeus 1758) based on capture-recapture data in a Venice Lagoon wetland area, Italy. Ethology Ecology & Evolution, 33(6), 561–575. https://doi.org/10.1080/03949370.2020.1870567 Lohr, S. C., P. A. Byorth, C. M. Kaya, and W. P. Dwyer. 1996. High temperature tolerances of fluvial Arctic grayling and comparisons with summer water temperatures of the Big Hole River, Montana. Transactions of the American Fisheries Society 125:933-939. Macnaughton, C. J., Durhack, T. C., Mochnacz, N. J., & Enders, E. C. (2021). Metabolic performance and thermal preference of westslope cutthroat trout (Oncorhynchus clarkii lewisi) and non-native trout across an ecologically relevant range of temperatures 1. Canadian Journal of Fisheries and Aquatic Sciences, 78(9), 1247–1256. https://doi.org/10.1139/cjfas-2020-0173 Martin, T. L., and R. B. Huey. 2008. Why “Suboptimal’’ is optimal: Jensen’s inequality and ectotherm thermal preferences. American Naturalist 171:E102–E118. Martins, E. G., Hinch, S. G., Patterson, D. A., Hague, M. J., Cooke, S. J., Miller, K. M., Lapointe, M. F., English, K. K., & Farrell, A. P. (2011). Effects of river temperature and climate warming on stock-specific survival of adult migrating Fraser River sockeye salmon (Oncorhynchus nerka). Global Change Biology, 17(1), 99–114. https://doi.org/10.1111/j.1365-2486.2010.02241.x Martins, E., O’Connor, B., Bottoms, J., Clevenger, I., Smith, D., Auger-Methe, M., Power, M., Patterson, D., Shrimpton, M., Cooke, S., 2022. Spatial Ecology of Arctic Grayling in the Parsnip Core Area. Report prepared for the Fish and Wildlife Compensation Program – Peace Region, Prince George, BC. Project Number: PEA-F22-F-3388. Marques, T. A., Thomas, L., Martin, S. W., Mellinger, D. K., Jarvis, S., Morrissey, R. P., Ciminello, C.-A., & DiMarzio, N. (2012). Spatially explicit capture–recapture methods to estimate minke whale density from data collected at bottom-mounted hydrophones. Journal of Ornithology, 152(S2), 445–455. https://doi.org/10.1007/s10336-010-0535-7 Masser, M. P., & Neill, W. H. (1986). Routes of heat transfer in two teleosts, Ictalurus punctatus and Lepomis macrochirus. Environmental Biology of Fishes, 16(4), 321–324. https://doi.org/10.1007/BF00842988 McCullough, D. A., Bartholow, J. M., Jager, H. I., Beschta, R. L., Cheslak, E. F., Deas, M. L., Ebersole, J. L., Foott, J. S., Johnson, S. L., Marine, K. R., Mesa, M. G., Petersen, J. H., Souchon, Y., Tiffan, K. F., & Wurtsbaugh, W. A. (2009). Research in thermal biology: Burning questions for coldwater stream fishes. Reviews in Fisheries Science, 17(1), 90–115. https://doi.org/10.1080/10641260802590152 McGarvey, R., Linnane, A. J., Feenstra, J. E., Punt, A. E., & Matthews, J. M. (2010). Integrating recapture-conditioned movement estimation into spatial stock assessment: A South Australian 176 lobster fishery application. Fisheries https://doi.org/10.1016/j.fishres.2010.03.006 Research, 105(2), 80–90. McPhail, J. D., & McPhail, D. L. (2007). The freshwater fishes of British Columbia. University of Alberta Press. Melnychuk, M. C., Dunton, K. J., Jordaan, A., McKown, K. A., & Frisk, M. G. (2017). Informing conservation strategies for the endangered Atlantic sturgeon using AT and multi-state markrecapture models. Journal of Applied Ecology, 54(3), 914–925. https://doi.org/10.1111/13652664.12799 Morash, A. J., Speers‐Roesch, B., Andrew, S., & Currie, S. (2021). The physiological ups and downs of thermal variability in temperate freshwater ecosystems. Journal of Fish Biology, 98(6), 1524–1535. https://doi.org/10.1111/jfb.14655 Mortensen, A., Ugedal, O., & Lund, F. (2007). Seasonal variation in the temperature preference of Arctic charr (Salvelinus alpinus). Journal of Thermal Biology, 32(6), 314–320. https://doi.org/10.1016/j.jtherbio.2007.03.004 Murphy, S. M., Adams, J. R., Waits, L. P., & Cox, J. J. (2021). Evaluating otter reintroduction outcomes using genetic spatial capture–recapture modified for dendritic networks. Ecology and Evolution, 11(21), 15047–15061. https://doi.org/10.1002/ece3.8187 Nakano, S., Kitano, S., Nakai, K., & Fausch, K. D. (1998). Competitive interactions for foraging microhabitat among introduced brook charr, Salvelinus fontinalis, and native bull charr, Bull trout Salvelinus confluentus, and westslope cutthroat trout, Oncorhynchus clarki lewisi, in a Montana stream. Environmental Biology of Fishes, 52(1–3), 345–355. https://doi.org/10.1023/A:1007359826470 Naman, S. M., Rosenfeld, J. S., Third, L. C., & Richardson, J. S. (2017). Habitat-specific production of aquatic and terrestrial invertebrate drift in small forest streams: Implications for drift-feeding fish. Canadian Journal of Fisheries and Aquatic Sciences, 74(8), 1208–1217. https://doi.org/10.1139/cjfas-2016-0406 Niella Y, Flávio H, Smoothey A, Aarestrup K, Taylor M, Peddemors V, Harcourt R (2020). “Refined Shortest Paths (RSP): incorporation of topography in space use estimation from node‐ based telemetry data.” Methods in Ecology and Evolution_. . O’Connor, B. 2023. Linking spatial stream network modeling and telemetry data to investigate thermal habitat use by adult arctic grayling. doi:https://doi.org/10.24124/2023/59378 Parlin, A. F., & Schaeffer, P. J. (2022). Cardiovascular contributions and energetic costs of thermoregulation in ectothermic vertebrates. Journal of Experimental Biology, 225(Suppl_1), jeb243095. https://doi.org/10.1242/jeb.243095 Pebesma, E., 2018. Simple Features for R: Standardized Support for Spatial Vector Data. The R Journal 10 (1), 439-446, https://doi.org/10.32614/RJ-2018-009 Pebesma, E., & Bivand, R. (2023). Spatial Data Science: With Applications in R (1st ed.). Chapman and Hall/CRC. https://doi.org/10.1201/9780429459016 177 Penney, C. M., Tabh, J. K. R., Wilson, C. C., & Burness, G. (2022). Within-generation and transgenerational plasticity of a temperate salmonid in response to thermal acclimation and acute temperature stress. Physiological and Biochemical Zoology, 95(6), 484–499. https://doi.org/10.1086/721478 Pépino, M., Goyer, K., Magnan, P. (2015); Heat transfer in fish: are short excursions between habitats a thermoregulatory behaviour to exploit resources in an unfavourable thermal environment?. J Exp Biol 1 November 2015; 218 (21): 3461–3467. doi: https://doiorg.prxy.lib.unbc.ca/10.1242/jeb.126466 Pincebourde, S., Murdock, C. C., Vickers, M., & Sears, M. W. (2016). Fine-scale microclimatic variation can shape the responses of organisms to global change in both natural and urban environments. Integrative and Comparative Biology, 56(1), 45–61. https://doi.org/10.1093/icb/icw016 Pirotta, E., Thompson, P. M., Cheney, B., Donovan, C. R., & Lusseau, D. (2015). Estimating spatial, temporal and individual variability in dolphin cumulative exposure to boat traffic using spatially explicit capture–recapture methods. Animal Conservation, 18(1), 20–31. https://doi.org/10.1111/acv.12132 Poole, G. C. (2002). Fluvial landscape ecology: Addressing uniqueness within the river discontinuum. Freshwater Biology, 47(4), 641–660. https://doi.org/10.1046/j.13652427.2002.00922.x Pörtner, H. O. (2002). Climate variations and the physiological basis of temperature dependent biogeography: Systemic to molecular hierarchy of thermal tolerance in animals. Comparative Biochemistry and Physiology Part A: Molecular & Integrative Physiology, 132(4), 739–761. https://doi.org/10.1016/S1095-6433(02)00045-4 Pörtner, H. O., & Farrell, A. P. (2008). Physiology and climate change. Science, 322(5902), 690– 692. https://doi.org/10.1126/science.1163156 Qasem, L., Cardew, A., Wilson, A., Griffiths, I., Halsey, L. G., Shepard, E. L. C., Gleiss, A. C., & Wilson, R. (2012). Tri-axial dynamic acceleration as a proxy for animal energy expenditure; should we be summing values or calculating the vector? PLoS ONE, 7(2), e31187. https://doi.org/10.1371/journal.pone.0031187 Raabe, J. K., Gardner, B., & Hightower, J. E. (2014). A spatial capture–recapture model to estimate fish survival and location from linear continuous monitoring arrays. Canadian Journal of Fisheries and Aquatic Sciences, 71(1), 120–130. https://doi.org/10.1139/cjfas-2013-0198 Rantanen, M., Karpechko, A.Y., Lipponen, A. et al. The Arctic has warmed nearly four times faster than the globe since 1979. Commun Earth Environ 3, 168 (2022). https://doi.org/10.1038/s43247022-00498-3 Raby, G. D., Donaldson, M. R., Hinch, S. G., Patterson, D. A., Lotto, A. G., Robichaud, D., English, K. K., Willmore, W. G., Farrell, A. P., Davis, M. W., & Cooke, S. J. (2012). Validation of reflex indicators for measuring vitality and predicting the delayed mortality of wild coho salmon bycatch released from fishing gears. Journal of Applied Ecology, 49(1), 90–98. https://doi.org/10.1111/j.1365-2664.2011.02073.x 178 Reid, A. J., Carlson, A. K., Creed, I. F., Eliason, E. J., Gell, P. A., Johnson, P. T. J., Kidd, K. A., MacCormack, T. J., Olden, J. D., Ormerod, S. J., Smol, J. P., Taylor, W. W., Tockner, K., Vermaire, J. C., Dudgeon, D., & Cooke, S. J. (2019). Emerging threats and persistent conservation challenges for freshwater biodiversity. Biological Reviews, 94(3), 849–873. https://doi.org/10.1111/brv.12480 Reid, C.H., C.S. Vandergoot, J.D. Midwood, E.D. Stevens, J. Bowker and S.J. Cooke. 2019. On the electroimmobilization of fishes for research and practice: opportunities, challenges, and research needs. Fisheries. 44(12):576-585. Rezende, E. L., Castañeda, L. E., & Santos, M. (2014). Tolerance landscapes in thermal ecology. Functional Ecology, 28(4), 799–809. https://doi.org/10.1111/1365-2435.12268 Riato, L., Leibowitz, S. G., Weber, M. H., & Hill, R. A. (2023). A multiscale landscape approach for prioritizing river and stream protection and restoration actions. Ecosphere, 14(1), e4350. https://doi.org/10.1002/ecs2.4350 Ritter, T. D., Zale, A. V., Grisak, G., & Lance, M. J. (2020). Groundwater upwelling regulates thermal hydrodynamics and salmonid movements during high-temperature events at a montane tributary confluence. Transactions of the American Fisheries Society, 149, 600–619. https://doi.org/10.1002/tafs.10259 Rodgers, E. M., & Gomez Isaza, D. F. (2022). Stress history affects heat tolerance in an aquatic ectotherm (Chinook salmon, Oncorhynchus tshawytscha). Journal of Thermal Biology, 106, 103252. https://doi.org/10.1016/j.jtherbio.2022.103252 Rose, D. C., Sutherland, W. J., Amano, T., González‐Varo, J. P., Robertson, R. J., Simmons, B. I., Wauchope, H. S., Kovacs, E., Durán, A. P., Vadrot, A. B. M., Wu, W., Dias, M. P., Di Fonzo, M. M. I., Ivory, S., Norris, L., Nunes, M. H., Nyumba, T. O., Steiner, N., Vickery, J., & Mukherjee, N. (2018). The major barriers to evidence‐informed conservation policy and possible solutions. Conservation Letters, 11(5), e12564. https://doi.org/10.1111/conl.12564 Rowe, D. K., & Chisnall, B. L. (1995). Effects of oxygen, temperature and light gradients on the vertical distribution of rainbow trout, Oncorhynchus mykiss, in two North Island, New Zealand, lakes differing in trophic status. New Zealand Journal of Marine and Freshwater Research, 29(3), 421–434. https://doi.org/10.1080/00288330.1995.9516676 Royle, J. A., Chandler, R. B., Sollmann, R., & Gardner, B. (Eds.). (2014). Spatial capturerecapture. Elsevier. Samuel, W. T., Yancy, L. E., Hinkle, E. G., & Falke, J. A. (2024). Validating morphometrics as a nonlethal tool to determine Arctic Grayling sex. North American Journal of Fisheries Management, 44(1), 70–78. https://doi.org/10.1002/nafm.10956 Schneider, D. C. (2001). The rise of the concept of scale in ecology. BioScience, 51(7), 545. https://doi.org/10.1641/0006-3568(2001)051[0545:TROTCO]2.0.CO;2 Sears, M. W., & Angilletta, M. J. (2015). Costs and benefits of thermoregulation revisited: Both the heterogeneity and spatial structure of temperature drive energetic costs. The American Naturalist, 185(4), E94–E102. https://doi.org/10.1086/680008 179 Sears, M. W., Angilletta, M. J., Schuler, M. S., Borchert, J., Dilliplane, K. F., Stegman, M., Rusch, T. W., & Mitchell, W. A. (2016). Configuration of the thermal landscape determines thermoregulatory performance of ectotherms. Proceedings of the National Academy of Sciences, 113(38), 10595–10600. https://doi.org/10.1073/pnas.1604824113 Sears, M. W., Riddell, E. A., Rusch, T. W., & Angilletta, M. J. (2019). The World Still Is Not Flat: Lessons Learned from Organismal Interactions with Environmental Heterogeneity in Terrestrial Environments. Integrative and Comparative Biology, 59(4), 1049–1058. https://doi.org/10.1093/icb/icz130 Schielzeth, H. 2010. Simple means to improve the interpretability of regression coefficients: Interpretation of regression coefficients. Methods in Ecology and Evolution 1:103–113. Schlosser, I. J. (1991). Stream fish ecology: A landscape perspective. BioScience, 41(10), 704– 712. https://doi.org/10.2307/1311765 Shepard, E. L. C., Wilson, R. P., Rees, W. G., Grundy, E., Lambertucci, S. A., & Vosper, S. B. (2013). Energy landscapes shape animal movement ecology. The American Naturalist, 182(3), 298–312. https://doi.org/10.1086/671257 Stamford, M., Hagen, J., & Williamson, S. (2017). FWCP Arctic grayling Synthesis Report: Limiting factors, enhancement potential, and critical habitats for Arctic grayling in the Williston Reservoir watershed, and information gaps limiting potential conservation and enhancement actions. Stanford, J. A., Lorang, M. S., & Hauer, F. R. (2005). The shifting habitat mosaic of river ecosystems. SIL Proceedings, 1922-2010, 29(1), 123–136. https://doi.org/10.1080/03680770.2005.11901979 Stevens, E. D., & Fry, F. E. J. (1974). Heat transfer and body temperatures in non thermoregulatory teleosts. Canadian Journal of Zoology, 52(9), 1137–1143. https://doi.org/10.1139/z74-152 Stevenson, R. D. (1985). Body size and limits to the daily range of body temperature in terrestrial ectotherms. The American Naturalist, 125(1), 102–117. https://www.jstor.org/stable/2461610 Stewart, D.B., Mochnacz, N., Reist, J.D., and Carmichael, T.J. 2007. Fish life history and habitat use in the Northwest Territories: Arctic Grayling (Thymallus arcticus). Can. Manuscr. Rep. Fish. Aquat. Sci. 2797: 64 Su, C. H., and Ryu, D. (2015). Multi-scale analysis of bias correction of soil moisture. Hydrol. Earth Syst. Sci. 19, 17–31. doi: 10.5194/hess-19-17-2015 Sutherland, C., Fuller, A. K., & Royle, J. A. (2015). Modelling non‐Euclidean movement and landscape connectivity in highly structured ecological networks. Methods in Ecology and Evolution, 6(2), 169–177. https://doi.org/10.1111/2041-210X.12316 Sutton, R. J., Deas, M. L., Tanaka, S. K., Soto, T., & Corum, R. A. (2007). Salmonid observations at a Klamath River thermal refuge under various hydrological and meteorological conditions. River Research and Applications, 23(7), 775–785. https://doi.org/10.1002/rra.1026 180 Thompson, P.D., Vasquez, E.A., Gowing, I., Edgar, T., Neville, A. and Jones, A. (2021), Unmanned Aerial Vehicle Technology Proves an Effective and Efficient Technique for Identifying Critical Native Fish Habitat. North Am J Fish Manage. https://doi.org/10.1002/nafm.10567 Tiffan, K. F., Kock, T. J., Connor, W. P., Steinhorst, R. K., & Rondorf, D. W. (2009). Behavioural thermoregulation by subyearling fall (Autumn) Chinook salmon Oncorhynchus tshawytscha in a reservoir. Journal of Fish Biology, 74(7), 1562–1579. https://doi.org/10.1111/j.10958649.2009.02228.x Tingley, R.W., Infante, D.M., Dean, E.M. et al. A landscape approach for identifying potential reestablishment sites for extirpated stream fishes: an example with Arctic grayling (Thymallus arcticus) in Michigan. Hydrobiologia 849, 1397–1415 (2022). https://doiorg.prxy.lib.unbc.ca/10.1007/s10750-021-04791-8 Torgersen, C. E., Price, D. M., Li, H. W., & McIntosh, B. A. (1999). Multiscale thermal refugia and stream habitat associations of chinook salmon in northeastern oregon. Ecological Applications, 9(1), 301–319. https://doi.org/10.1890/10510761(1999)009[0301:MTRASH]2.0.CO;2 Torgersen, C. E., Le Pichon, C., Fullerton, A. H., Dugdale, S. J., Duda, J. J., Giovannini, F., Tales, É., Belliard, J., Branco, P., Bergeron, N. E., Roy, M. L., Tonolla, D., Lamouroux, N., Capra, H., & Baxter, C. V. (2022). Riverscape approaches in practice: Perspectives and applications. Biological Reviews, 97(2), 481–504. https://doi.org/10.1111/brv.12810 Valle, D., Mintz, J., & Brack, I. V. (2024). Estimation and interpretation problems and solutions when using proportion covariates in linear regression models. Ecology, 105(4), e4256. https://doi.org/10.1002/ecy.4256 Vannote, R. L., Minshall, G. W., Cummins, K. W., Sedell, J. R., & Cushing, C. E. (1980). The river continuum concept. Canadian Journal of Fisheries and Aquatic Sciences, 37(1), 130–137. https://doi.org/10.1139/f80-017 Venables, W. N. & Ripley, B. D. (2002) Modern Applied Statistics with S. Fourth Edition. Springer, New York. ISBN 0-387-95457-0 Vickers, M., Manicom, C., Schwarzkopf, L., Demas, A. E. G. E., & McPeek, E. M. A. (2011). Extending the cost-benefit model of thermoregulation: High-temperature environments. The American Naturalist, 177(4), 452–461. https://doi.org/10.1086/658150 Veech, J. A. (2021). Habitat ecology and analysis. Oxford university press. Wagner, G. F., S. J. Cooke, R. S. Brown, and K. A. Deters. 2011. Surgical implantation techniques for electronic tags in fish. Reviews in Fish Biology and Fisheries 21:71–81. Webb, J. K., & Whiting, M. J. (2005). Why don’t small snakes bask? Juvenile broad‐headed snakes trade thermal benefits for safety. Oikos, 110(3), 515–522. https://doi.org/10.1111/j.00301299.2005.13722.x Weins, J.A. (1989). Spatial scaling in ecology. Functional Ecology 3: 385–397. 181 Wellemeyer, J. C., Perkin, J. S., Jameson, M. L., Costigan, K. H., & Waters, R. (2019). Hierarchy theory reveals multiscale predictors of Arkansas darter (Etheostoma cragini) abundance in a Great Plains riverscape. Freshwater Biology, 64(4), 659–670. https://doi.org/10.1111/fwb.13252 Weller, D. E., Anderson, D. J., DeAngelis, D. L., & Coutant, C. C. (1984). Rates of heat exchange in largemouth bass: Experiment and model. Physiological Zoology, 57(4), 413–427. https://doi.org/10.1086/physzool.57.4.30163343 White, S. L., Kline, B. C., Hitt, N. P., & Wagner, T. (2019). Individual behaviour and resource use of thermally stressed brook trout Salvelinus fontinalis portend the conservation potential of thermal refugia. Journal of Fish Biology, 95(4), 1061–1071. https://doi.org/10.1111/jfb.14099 Whoriskey, K., Martins, E. G., Auger‐Méthé, M., Gutowsky, L. F. G., Lennox, R. J., Cooke, S. J., Power, M., & Mills Flemming, J. (2019). Current and emerging statistical techniques for aquatic telemetry data: A guide to analysing spatially discrete animal detections. Methods in Ecology and Evolution, 10(7), 935–948. https://doi.org/10.1111/2041-210X.13188 Wich, S. A., & Koh, L. P. (2018). Conservation drones (Vol. 1). Oxford University Press. https://doi.org/10.1093/oso/9780198787617.001.0001 Williams J.E. 2000. The Coefficient of Condition of Fish. Chapter 13 in Schneider, James C. (ed.) 2000. Manual of fisheries survey methods II: with periodic updates. Michigan Department of Natural Resources, Fisheries Special Report 25, Ann Arbor. Wilms, T., & Whitworth, G. (2016). Mapping of critical summer thermal refuge habitats for chinook salmon, coho salmon, steelhead and bull trout in the Nicola River Watershed – 2016 (Habitat Stewardship Program for Species at Risk project reference number: 2016HSP7592). Wojcik, F. J. 1955. Life history and management of the grayling in interior Alaska. Master’s thesis. University of Alaska. Wolkovich, E. M., Cook, B. I., McLauchlan, K. K., & Davies, T. J. (2014). Temporal ecology in the Anthropocene. Ecology Letters, 17(11), 1365–1379. https://doi.org/10.1111/ele.12353 Woodget, A. S., Austrums, R., Maddock, I. P., & Habit, E. (2017). Drones and digital photogrammetry: From classifications to continuums for monitoring river habitat and hydromorphology. WIREs Water, 4(4), e1222. https://doi.org/10.1002/wat2.1222 Zemlak R. J., and A. R. Langston. 1998. Fish species presence and abundance of the Table River, 1995. Peace/Williston Fish and Wildlife Compensation Program PWFWCP Report No. 173. Zhao, B., Ding, P., & Zhang, W. (2014). Energetic cost of behavioral thermoregulation in turtle embryos. Behavioral Ecology, 25(4), 924–927. https://doi.org/10.1093/beheco/aru077 Zuur, A. F., Ieno, E. N., Walker, N., Saveliev, A. A., & Smith, G. M. (2009). Mixed effects models and extensions in ecology with R. Springer New York. https://doi.org/10.1007/978-0-387-874586 182 Supplemental Information Summer: D ~ Temp + Temp2 + BT Model parameters Estimate SE D (intercept) Temp Temp2 BT α σ Model parameters D (intercept) BT α σ Model parameters D (intercept) α σ -2.92 -1.12 -0.57 0.32 -2.50 9.07 Spring: D ~ BT Estimate -2.75 0.18 -2.89 8.75 Winter: D ~ 1 Estimate -2.23 -1.65 7.11 lcl ucl 0.19 0.32 0.23 0.10 0.06 0.03 -3.30 -1.76 -1.02 0.13 -2.61 9.02 -2.54 -0.49 -0.12 0.51 -2.39 9.12 SE 0.14 0.08 0.09 0.04 lcl -3.02 0.02 -3.06 8.66 ucl -2.48 0.34 -2.71 8.83 SE 0.13 0.07 0.02 lcl -2.48 -1.78 7.07 ucl -1.98 -1.51 7.16 SI 2.1. Model coefficients for the top T. arcticus models for each season. Covariates Temp, Temp2, and BT represent the polynomial temperature coefficients and the predator S. confluentus covariate, α is the detection probability, and σ is the scale parameter. The standard error and upper and lower confidence intervals are given as SE, ucl, and lcl, respectively. 183 α σ D σ α D σ α D 0.00 2021 Estimate SE -2.90 0.22 -3.53 0.06 9.51 0.04 18.99 0.06 0.04 2020 Estimate SE -2.91 0.22 -6.25 0.06 -3.53 9.51 18.99 ucl -2.46 -3.41 9.58 lcl -3.34 -3.66 9.44 ucl -2.47 -6.13 lcl -3.35 -6.37 18.99 -3.41 9.58 ucl -2.46 -3.66 9.44 Summer 2019 Estimate SE lcl -2.90 0.22 -3.34 α σ D σ α D σ α D 0.09 2021 Estimate SE -4.73 0.58 -0.97 0.25 8.29 0.08 8.55 0.64 0.38 2020 Estimate SE -4.04 0.41 -2.08 0.20 -3.82 8.84 Spring 2019 Estimate SE -5.01 0.72 lcl -5.86 -1.47 8.13 8.38 lcl -4.84 -2.47 -5.08 8.10 lcl -6.42 ucl -3.60 -0.48 8.46 8.72 ucl -3.24 -1.69 -2.56 9.59 ucl -3.60 α σ D σ α D σ α D 0.05 2021 Estimate SE -4.16 0.45 -2.81 0.20 7.69 0.09 8.46 0.20 0.06 2020 Estimate SE -3.53 0.32 -3.17 0.11 -1.72 7.63 Winter 2019 Estimate SE -4.69 0.58 lcl -5.04 -3.21 7.52 8.35 lcl -4.15 -3.39 -2.11 7.52 lcl -5.83 184 ucl -3.29 -2.41 7.87 8.56 ucl -2.91 -2.95 -1.34 7.75 ucl -3.56 SI 2.2. Model coefficients for S. confluentus models. There were no covariate effects in this model, so the response parameters are the intercept, detection probability α, and the scale parameter σ. The standard error and upper and lower confidence intervals are given as SE, ucl, and lcl, respectively. 185 SI 2.3. Seasonal distance of tagged T. arcticus activity centres from their location of tagging. Tags are disproportionately near their tagging sites in the summer season, as T. arcticus show high stream fidelity and they were tagged in the summer, though some individuals migrated further afield (up to 73 Rkm away). Spring and winter activity centres were predicted relatively evenly between zero and 52 Rkm away from the tagging sites. 186 SI 2.4. The mean fork length in millimetres of S. confluentus detected by the acoustic receiver array during each season. Smaller fish were detected during the winter season, likely a result of a component large adfluvial fish within the population which leaves the study area to overwinter in the Williston Reservoir. SI 2. 5. Seasonal distributions of the S. confluentus covariate BT across the Parsnip watershed. 187 188 SI 2.6. Acoustic receivers which were lost or damaged due to high water flows. Particularly notable losses occurred in the Anzac and Table Rivers during the high discharge freshet and sustained high flows in the 2019 and 2020 seasons. Acoustic receivers were replaced at suitable sites when water levels returned to working conditions. 189 SI 3.1. Arctic grayling counts at both the reach scale (2022, 100-m resolution) and river scale (2018, 2019, 2021 at a 1.2 – 5.4 Rkm resolution). Data for the reachscale surveys were produced for this study and data for the riverscape study were provided by John Hagen and Associates (Hagen and Stamford 2023). Date 2022-08-11 2022-08-11 2022-08-11 2022-08-12 2022-08-12 2022-08-12 2022-08-12 2022-08-14 2022-08-14 2022-08-14 2022-08-14 2022-08-14 2022-08-14 Flight 1 2 3 4 5 6 7 8 9 10 11 12 13 15:53 14:40 13:43 12:49 11:59 10:49 16:48 16:04 14:33 10:38 15:07 13:38 13:09 Start 16:09 14:55 13:57 12:58 12:13 11:02 17:09 16:12 14:45 10:53 15:25 13:50 13:21 Stop SI 3.2. Summary information for the river-scale drone flights. 190 191 SI 3.3. Example of a composite image built for this study which includes the RGB raster (panel A), the FLIR raster (Panel B), and a 0.25 °C contour layer (Panel C) which clarifies the thermal distributions within the imagery. White crosses and labels represent point sample temperature measurements from the FLIR layer. 192 The workflows used in this study to produce and analyze drone imagery were produced by manual-flight mapping in lieu of preprogrammed flights due to its relative efficiency compared with much slower pre-programmed flights when working at the spatial scales I explore in this chapter. This approach comes with tradeoffs, as image orthorectification can be sensitive to variable flight heights, speeds, and image overlap (Wich and Pin Koh 2018; Dugdale 2016; Dugdale 2015). Using Pix4DMapper’s photogrammetry tools which accommodate variable image overlap, we were successful in producing high-quality orthomosaics at both the reach and river scales at flight speeds up to 32 kph and heights up to 100 m). As FLIR sensors such as the one used in this study only measure the surface temperatures of an object, it is possible that deeper pools had more internal thermal heterogeneity than what was represented by surface signatures, though the continuous mixing of rivers assumes that surface temperatures are largely representative of the water column beneath (Dugdale 2016). To extract only wetted pixels for use in the reach-scale analyses, river masks were created by rendering contour lines from the FLIR layer and assembling a polygonised mask layer from isotherms selected at the land-water interface. Minor spatial differences in the completed RGB and FLIR rasters (derived from the different pixel widths and resolutions of the dual cameras aboard the imaging payload) were corrected by georeferencing the RGB layers to the FLIR layers using a spline transformation with at least 10 control points per layer in ArcGIS Pro (version 3.1.2; Advanced License; ESRI, Redlands, CA, U.S.A.). By georeferencing the layers in this order (i.e. RGB to FLIR), the spatial structure of the radiometrically sensed pixels of the FLIR raster was preserved. RGB and FLIR rasters were then imported into QGIS software (version 3.32.1; GNU General Public License; London, England) for further processing to create 0.25 °C contours and assemble composite image layers for visual analysis (Figure 5). In Pix4Dmapper by Pix4D photogrammetry software (version 4.8.3; Educational License; Pix4D S.A., Prilly, Switzerland), image pair matching was set to Free Flight or Terrestrial, and FLIR rasters and RGB orthomosaics export options were set to merged GeoTIFF files. All other point cloud processing, camera configuration, and image georeferencing and datum selection options were left as Automatic and were derived from EXIF data stored in the raw drone image files. SI 3.4. Expanded description of drone processing methods. 193 SI 3.5. Covariate spatial summaries used in the river-scale models. The top three panels represent the SSNM temperatures (O'Connor 2023) for 2018, 2019, and 2021. The bottom two panels represent mean pool size and pool densities at a 1 Rkm aggregation. SI 3.6. Habitat classes defined for the reach-scale study. Polygons were defined manually by tracing habitat features over RGB drone orthomosaics. 194 195 SI 3.7. Model averaged predictions of the effects of bull trout counts (Panel A) and mountain whitefish presence or absence (Panel B) on Arctic grayling abundance at the reach scale. Arctic grayling abundance decreased with bull trout abundance and increased with mountain whitefish presence. -5.44 0.05 1.77 2.35 hab.pool x temp hab.pool x temp2 Dispersion Standard deviation of random effect -2.41 1.28 temp temp 0.36 hab.pool 2 0.45 -0.10 0.64 Estimate (Intercept) BT MW Term - - 1.41 1.34 1.40 1.28 0.41 0.46 0.16 0.49 Std. Error - - 0.03 -4.06 -1.72 1.00 0.89 0.96 -0.61 1.32 z value - - 0.97 0.00 0.08 0.32 0.37 0.34 0.54 0.19 Pr(>|z|) 196 SI 3.8. Summary table for reach-scale global model coefficients. Terms are bull trout (b), mountain whitefish (w), area of pool habitats in a subreach (h.pool), and the polynomial temperature term (t + t 2). Bold terms indicate which of the two polynomial temperature terms the model corresponds to. 197 SI 3.9. Quantile residuals for reach-scale global model. The left plot shows the QQ plot of residuals versus expected values, and the right plot shows the spread of the residuals by quantile. 198 SI 3.10. Visual assessment of residuals by Rkm to assess for spatial autocorrelation from reach-scale global model. The vertical line is a partition for the data from reach A (left of the line) and reach B (right of the line). 199 SI 3.11. RGB imagery from the four survey dates in reach A. Wetted width remained approximately the same between sampling dates (relative to freshet and inseason discharge spikes observed during other years in this system), though water levels generally dropped throughout the sampling season. Beginning in the August 15 imagery and becoming most pronounced on September 2 baseflows before autumn rainfalls began, gravel bars begin emerging in the middle and lower portions of the reach, the lower of which persisted even after rainfall boosted water levels in the late season. 200 SI 3.12. RGB imagery from the four survey dates in reach B. In comparison to reach A, channel width constrictions in reach B due to dropping water levels was slightly more apparent, but still not largely impacted. Vertical water levels impacted clarity more prominently between images in reach B than in reach A. As water levels dropped in reach B, an ephemeral side channel in the middle of the reach dried up and did not return after the July 31 survey. 7 0 0 0 GR BT MW GR 0 0 MW MW 0 BT 1 7 GR AR02 0 0 13 0 0 3 1 1 1 1 AR03 BT 0 MW 1 0 2 1 5 7 1 0 0 1 AR04 0 0 4 1 2 4 1 1 2 1 1 AR05 1 0 0 0 0 0 0 0 0 0 0 0 AR06 2 0 0 15 0 0 4 0 1 1 0 0 AR07 A4 (Sep. 15) 2 BT 1 1 2 0 0 3 1 4 1 1 1 AR08 1 0 0 1 0 0 1 3 2 1 0 0 0 0 4 0 0 5 0 0 0 0 0 0 AR09 5 - - - - - - - - - 0 0 0 AR10 0 - - - - - - - - - 0 0 0 AR11 14 - - - - - - - - - 0 0 0 AR12 1 - - - - - - - - - - - - AR13 3 0 0 0 0 0 1 0 0 0 0 0 0 AR14 18 0 0 0 0 0 0 0 0 1 0 0 2 AR15 A3 (Sep. 2) 7 GR 0 0 0 - - - 0 0 0 0 0 0 AR16 A2 (Aug. 15) AR01 Spp. 0 0 0 0 0 0 0 0 0 0 3 0 AR17 A1 (Aug. 2) Survey Reach A Snorkel Survey Counts by Species 3 2 40 3 7 34 5 10 15 5 10 50 201 SI 3.13. Counts of Arctic grayling (GR), bull trout (BT), and presence/absence of mountain whitefish (MW) from the reach-scale snorkel surveys in reach A. Reaches which could not be snorkeled due to water conditions are listed as no data (-). TOTAL 0 2 MW GR 0 0 BT MW 6 GR 0 0 MW BR02 0 0 0 0 0 0 - - BR03 BT 0 BT 1 0 2 1 2 10 1 0 BR04 0 0 0 0 0 0 0 0 0 0 6 0 0 0 0 0 32 BR05 9 - - - - - - - - - BR06 47 - - - - - - - - - BR07 - - - - - - - - - - BR08 B4 (Sep. 15) 3 GR 0 0 4 1 0 15 0 0 44 BR09 B3 (Sep. 2) BR01 Spp. 0 0 0 0 0 0 0 0 3 BTRIB B2 (Aug. 15) Survey Reach B Snorkel Survey Counts by Species 1 0 14 2 2 31 1 0 138 202 SI 3.14. Counts of Arctic grayling (GR), bull trout (BT), and presence/absence of mountain whitefish (MW) from the reach-scale snorkel surveys in reach B. Reaches which could not be snorkeled due to water conditions are listed as no data (-). TOTAL 3.44 0.21 0.09 0.53 13.5 Pool density Rkm Dispersion Estimate (Intercept) Mean pool size Term SI 3.15. Summary table for river-scale global model coefficients. 0.07 0.06 0.06 0.06 Std. Error 7.47 1.43 58.57 3.41 z value 0.00 0.15 0.00 0.00 Pr(>|z|) 203 204 SI 3.16. Quantile residuals for reach-scale global model. The left plot shows the QQ plot of residuals versus expected values, and the right plot shows the spread of the residuals by quantile. SI 3.17. Visual assessment of residuals by Rkm to assess for spatial autocorrelation from reach-scale global model. 205 SI 3.18. Summary plots of mean pool size (Panel A) and pool density (Panel B) at a 1 Rkm aggregation. Pool size was defined by the pool scoring criteria described in text and in Figure 7 and pool density was the total count of pools within each reach. 206 SI 3.19. Modeled temperatures by river kilometer (Rkm) for 2018, 2019, and 2021. 207 208 SI 4.1. Plot illustrating the simulation approach to determining body temperature with 100 draws from a normal distribution of with mean = ܶ஺ and standard deviation of 1.5 °C. The coefficient of variation (Tb_cv) was used to cut data to less than 0.05 at which points body temperatures converged, or at which point the variation visible in the plot stabilized. 209 SI 4.2. Still images from the GIFs showing the hourly interpolation of the drone rasters produced in chapter 2 used to calculate the metric of the patchiness of the thermal habitat used in the GAMMs. Animated versions are available online for reach A and reach B. Sex F F M F F M M F F M Experiment ID 3 4 5 6 7 8 9 10 11 12 355 345 304 343 293 362 330 332 353 312 Fork length (mm) 500 500 450 500 350 540 425 490 660 360 Weight (g) A A A B A B B A B B Capture reach 11.8 12.4 12.0 13.6 14.3 13.0 12.3 15.1 13.7 14.3 TACC (°C) 12.6 13.9 9.0 9.0 7.2 18.4 14.5 11.3 7.6 12.6 TPREF (°C) 210 10.28 - 15.2 9.18 - 17.03 7.98 - 9.86 7.4 - 9.67 6.16 - 9.59 17.55 - 18.82 13.14 - 15.5 10.5 - 12.84 6.97 - 8.39 11.59 - 13.24 TSET (°C) SI 4.3. Table of acclimation temperatures ܶ஺஼஼ , or the mean river temperatures at the capture location for 2 weeks prior to the start of the shuttlebox thermal preference experiments. There was weak negative correlation between ܶ஺஼஼ and ܶ௉ோாி (-0.33, P = 0.36). SI 4.4. Plot of the effects of centred body weight (Wc) and the trend (warm) factor. The cooling trend is represented by the intercept, and warming trend effect estimates are shown for the top 2 models by AICc weight (m1 and m2, respectively). 211 SI 4.5. A pair of synchronized Star ODDI nano DST temperature loggers prepared for internal and external application. 212 Time 13:56 14:09 14:20 14:30 14:42 14:45 15:09 13:37 13:46 14:09 14:30 14:46 15:02 15:16 15:45 18:19 18:37 18:47 18:56 13:44 Date 2022-07-28 2022-07-28 2022-07-28 2022-07-28 2022-07-28 2022-07-28 2022-07-28 2022-08-05 2022-08-05 2022-08-05 2022-08-05 2022-08-05 2022-08-05 2022-08-05 2022-08-05 2022-08-05 2022-08-05 2022-08-05 2022-08-05 2022-08-06 A B B B B B B B B B B B B A A A A A A A Reach 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 Fish no. F M M F F M F F F F M F M M F F F M M F Sex 309 316 324 310 352 350 346 345 340 321 395 318 380 340 314 345 359 365 361 330 Len. (mm) 340 400 400 375 575 450 450 450 490 300 700 350 650 490 320 420 430 430 420 415 Wt. (g) 13 18 6 2 4 1 28 26 16 7 15 17 8 25 27 20 24 21 22 19 Radio no. 213 SI 4.6. Summary information for the radio-tagged Arctic grayling included in the index E analyses. Italicized terms for sex represent fish which were estimated as male or female based on the ratio of dorsal length to fork length. The reach column denotes where the individual was tagged, fish no. denotes the sequential number across all fish handled in this study, and the radio no. denotes the radio tag ID of the fish. Rows in bold indicate fish that were included in the index E analyses as they had sufficient data to pass the filtering steps. 18:28 19:11 20:52 21:04 12:27 12:33 12:39 12:45 17:13 17:18 17:22 17:35 17:42 17:53 20:07 20:19 15:33 15:54 20:20 20:25 20:33 20:40 13:52 14:01 14:14 14:21 2022-08-06 2022-08-06 2022-08-06 2022-08-06 2022-08-07 2022-08-07 2022-08-07 2022-08-07 2022-08-13 2022-08-13 2022-08-13 2022-08-13 2022-08-13 2022-08-13 2022-08-13 2022-08-13 2022-08-13 2022-08-13 2022-08-14 2022-08-14 2022-08-14 2022-08-14 2022-08-21 2022-08-21 2022-08-21 2022-08-21 B B B B B B B B B B B B B B B B B B B B B B B B B A 46 45 44 43 42 41 40 39 38 37 36 35 34 33 32 31 30 29 28 27 26 25 24 23 22 21 F M M F F M M F F F F M F F M M F M M F M M F F M F 299 329 398 308 335 325 305 344 307 380 390 328 324 350 375 383 350 340 313 310 338 337 354 310 408 323 350 400 625 350 460 420 390 510 340 380 440 400 440 500 640 580 560 390 640 620 799 640 780 610 800 675 46 43 42 48 41 40 39 38 37 36 34 33 30 50 29 35 32 31 23 9 3 10 11 14 12 5 214 14:30 15:18 15:27 15:37 2022-08-21 2022-08-21 2022-08-21 2022-08-21 B B B B 50 49 48 47 F F M M 335 328 389 370 600 420 600 550 47 49 45 44 215 Appendix A. SCR Modeling tutorial This appendix has been formatted for hosting online as a GitBook. It cannot be published live until the manuscript has been accepted. Minor formatting issues in this version are attributed to this. 216 Linking acoustic telemetry data to spatial covariates in dendritic river networks with spatially explicit capture-recapture (SECR) models - Tutorial Joseph R. Bottoms, Marie Auger-Méthé, Bryce O’Connor, Michael Power, David A. Patterson, Mark Shrimpton, Steven J. Cooke, and Eduardo G. Martins 2023-12-29 Introduction The purpose of this document is to demonstrate how to link acoustic telemetry data to spatially-explicit covariates in a branching (dendritic) river network using spatially-explicit capture-recapture (SECR) models. This tutorial uses preloaded data objects, inspects the data structures necessary for the analysis, and uses the secr and secrlinear packages (Efford 2023a,b) to fit, select, predict, and visualize the models. Our case study uses acoustic telemetry detections of Arctic grayling (Thymallus arcticus) from 2019-2021 in the Parsnip watershed in north-central British Columbia, Canada. For demonstration purposes, we isolate only the summer seasons for the tutorial. This document accompanies the manuscript Bottoms et al. (in prep). Due to the sensitivity of presenting spatial data on the distributions of protected bull trout (Salvelinus conflulentus), this tutorial demonstrates how to fit models using the temperature covariate from our case study and does not include the data used to create the bull trout covariate. For tutorial purposes, one covariate will suffice to demonstrate the workflow. The secrlinear package (Efford 2023b) is a wrapper for the parent secr package (Efford 2023a) that enables SECR models to be fit in linear state-spaces such as dendritic river networks. Applying this approach to a 2D open water system would be similar, though care should be taken to understand how data preparation, model specifications, and computation times differ in 2D systems. Please refer to the vignettes for the secr package (Efford 2023a) for more. This workflow produces a prediction of how the temperature covariate influences the seasonal distribution of Arctic grayling in the Parsnip watershed. While not explicitly a model of resource selection, it is most intuitively understood as a prediction of how individuals use the temperature resource within their summer home range. This model assumes demographic closure, which implies that all tags detected in a given window remain alive at the end of the window and does not account for inseason mortality or tag loss. It further assumes that there is only one activity centre (average location) for each tag during each seasonal window; fine-scale movements or movement parameters between seasons are not modeled by this approach. When applying this workflow, we recommend pre-analysis using visualization techniques provided in the actel (Flávio and Baktoft 2021) and RSP (Niella et al. 2020) packages to ensure that the windows being modeled are ecologically-relevant to the fine-scale distributions of the focal species and meet these assumptions. 217 For applications using acoustic telemetry data, we found SCR to be advantageous as it (a) uses the full detection history of the tagged population to link activity centres (average detected locations) of the tagged animals to spatially-explicit environmental covariates (analogous to resource selection at the home range scale), (b) can accommodate discontinuous acoustic receiver arrays, which is beneficial in systems where acoustic receiver loss and/or redeployment between study years is characteristic, and (c) pairs well with existing descriptive approaches to acoustic telemetry analysis to offer a statistically robust tool using a priori information derived from these methods about a species’ movement patterns. Due to the nature of acoustic telemetry datasets in which only the tagged population is available for detection, we highlight several specific considerations and assumptions for using this approach: i.e. (a) activity centres are fixed within modeling windows by the underlying closed population model, so this method is best applied to discrete, ecologicallyrelevant windows in which fine-scale movements are not the focus, (b) inferences from acoustic telemetry data depict relative (not absolute) densities of only the tagged population, and (c) spatial tagging effort must be defined in the model to ensure that predictions are not merely an artifact of tagging effort across space and time. # R version 4.2.3 library(secrlinear) library(secr) # Version 4.6.5 pre-release from GitHub MurrayEfford/sec r # devtools::install_github("MurrayEfford/secr") library(tidyverse) library(beepr) library(ggpubr) library(sf) library(kableExtra) Read and inspect prepared data files Read traps files. Separate traps objects are made for each year as acoustic receiver sites have varying usage between years. Refer to ?secr::read.traps() for the arguments required to create a traps objects. These objects were prepared using the arguments detector = 'proximity' and binary.usage = TRUE. traps.su19 <- readRDS('data/traps.su19.rds') traps.su20 <- readRDS('data/traps.su20.rds') traps.su21 <- readRDS('data/traps.su21.rds') Inspect the data structure of a traps object’s usage. Each trap (numbered row) is either active (1) or inactive (0) on each day/occasion (numbered column). In our case study, there were 76 unique acoustic receivers used across all years, but only a portion of these were active during each modeling window. The binary matrix which corresponds to a 218 modeling session can be assigned to a traps object using secr::usage(traps) "2019-07-01 to 2019-09-14", "2019-07-01 to 2019-09 -14", "2019… ## $ ID 24302, 24302, 24302, 24302, 24302, 24362, 24369, 2 4361, 24353… ## $ occasion 1, 2, 5, 5, 5, 8, 8, 8, 8, 9, 10, 10, 11, 11, 11, 11, 11, 12,… ## $ x 567807.0, 567807.0, 563837.8, 559575.8, 557700.5, 579001.3, 5… ## $ y 6061318, 6061318, 6058268, 6054844, 6054550, 60668 41, 6066841… Read covariate data Covariates must be available for each pixel in the state-space. They can also be defined at the tag level or the acoustic receiver level, but those are not explored in this tutorial. Here we have annual mean summer temperature values (rows) for each of the 366 pixels in the Parsnip watershed state-space. The continuous temperature covariate was created using imputed temperature data in a spatial stream network model (O’Connor 2023). covs <- readRDS('data/covariates.rds') glimpse(covs) ## Rows: 366 ## Columns: 5 ## $ x 529904.3, 530234.2, 530674.6, 531439.9, 531836.0, 5 31876.1, 53… ## $ y 6072795, 6072772, 6072335, 6071980, 6072850, 607379 9, 6074595,… ## $ temps19 11.68722, 11.68722, 11.70152, 11.66809, 11.62223, 1 1.57152, 11… ## $ temps20 9.734608, 9.734803, 9.731571, 9.732794, 9.708983, 9 .659940, 9.… ## $ temps21 11.32968, 11.32968, 11.32185, 11.31386, 11.29521, 1 1.28682, 15… ggplot() + geom_sf(data = covs %>% st_as_sf(coords = c('x','y')), aes(color = t emps20)) + theme_bw() 220 Standardize covariates for analysis covs <- covs %>% mutate(ztemp19 = (temps19 - mean(temps19)) / sd(temps19), ztemp20 = (temps20 - mean(temps20)) / sd(temps20), ztemp21 = (temps21 - mean(temps21)) / sd(temps21), ztemp2.19 = ztemp19^2, ztemp2.20 = ztemp20^2, ztemp2.21 = ztemp21^2) Create masks Masks are defined for each year, as each will have unique covariate values. habitatmap <- ("data/pcr.shp") # Multiline shapefile of the Parsnip Co re Region spacing <- 1000 # Pixel size of 1,000 m mask2019 <- read.linearmask(file = habitatmap, spacing = spacing) mask2020 <- read.linearmask(file = habitatmap, spacing = spacing) mask2021 <- read.linearmask(file = habitatmap, spacing = spacing) Covariate columns can be assigned to mask objects using covariates(mask)[,'cov.name'] <- cov.name. Since the multi-session model fits pooled data across all sessions, the covariate names must be the same in each sessional mask. Also refer to secr::addCovariates(). 221 covariates(mask2019)[,'temp'] <- covs$ztemp19 covariates(mask2020)[,'temp'] <- covs$ztemp20 covariates(mask2021)[,'temp'] <- covs$ztemp21 covariates(mask2019)[,'temp2'] <- covs$ztemp2.19 covariates(mask2020)[,'temp2'] <- covs$ztemp2.20 covariates(mask2021)[,'temp2'] <- covs$ztemp2.21 Sampling effort The distribution of tagging effort can bias model predictions if not accounted for. While a uniform tagging program should be applied to evenly distribute tags within the acoustic receiver array, logistical realities of sampling dynamic aquatic species and systems can at times lead to an uneven distribution of tags. To provide the model a way to account for tagging effort in space and time and to check that any relationships found between the tags and their spatial covariates are not simply an artifact of tagging effort, a spatial effort layer is created. This layer is then added into the model as a covariate which is used to thin the likelihood (Borchers and Efford 2008) by the probability that a tag detected at a given site had been previously tagged there. The probability surfaces used in this example were created by kernel density estimation using qgis_run_algorithm('qgis:heatmapkerneldensityestimation') from the qgisprocess package (Dunnington et al. 2023) with the kernel radius set to ߪ and the weight field defined as the cumulative tags applied at each tagging pixel divided by the total tags applied to-date (inclusive of the window being modeled). From this, we subtracted the tagging efforts greater than two years old (an approximation of the ~2.2 battery life of the tags). A better metric for this object would factor in catch per unit effort (CPUE), which would refine the probability surface to include hours spent tagging at sites which were both productive and unproductive to tagging. This example accumulates effort seasonally, which in the case of Arctic grayling which show high fidelity to their summer (tagging season) habitats was reasonable. Applications of this approach should consider the ecology of a species of interest when defining an effort layer and modeling windows. Here we visualize the summer effort layers used in each season. effort <- readRDS('data/effort.rds') ggplot() + geom_sf(data = effort, aes(color = gr.su19)) + theme_bw() + coord_sf(datum = st_crs(32610)) 222 ggplot() + geom_sf(data = effort, aes(color = gr.su20)) + theme_bw() + coord_sf(datum = st_crs(32610)) 223 ggplot() + geom_sf(data = effort, aes(color = gr.su21)) + theme_bw() + coord_sf(datum = st_crs(32610)) Assign effort covariates (pd for pdot in Efford’s terminology) to each mask covariates(mask2019)$pd <- effort$gr.su19 covariates(mask2020)$pd <- effort$gr.su20 covariates(mask2021)$pd <- effort$gr.su21 Define capture history Define the detection history of Arctic grayling as a multisession capthist object. Use secr::verify to check the integrity of the data or highlight data which is in conflict. GR.su <- make.capthist(captures = det.dataGR, traps = list(traps.su19, traps.su20, traps.su21), fmt = 'XY', bysession = TRUE) verify(GR.su) ## No errors found :-) 224 Fit candidate SECR models to data Once all the data is loaded and prepared, SECR models can be fit using the arguments as defined by ?secr::secr.fit(). Model su.m0 is the null model, modeling D against the only tagging effort layer pd. Models su.m1 and su.m2 add terms for temp and its squared term temp2, respectively. Starting parameters: D is 1; g0 and sigma were defined after exploratory SECR modeling of acoustic telemetry detections using intercept-only models to inform appropriate starting values. Details include the use of networkdistance, which enables likelihood estimation in our dendritic river network and relativeD = TRUE, which is required when fitting models to acoustic telemetry data in which only the relative density of tagged animals can be estimated. The argument steptol is a tuning parameter for fitting the model and is not necessary in all cases. Call ?nlm() to view the maximization error codes produced by secr.fit() to get a sense of which tuning parameters to define if you run into convergence issues. su.m0 <- secr.fit(capthist = GR.su, mask = list(mask2019, mask2020, mask2021), model = D ~ pd, start = list(D = 1, g0 = 0.1, sigma = 9000), details = list(userdist = networkdistance, relativeD = TRUE), trace = FALSE, steptol = 1e-3) su.m1 <- secr.fit(capthist = GR.su, mask = list(mask2019, mask2020, mask2021), model = D ~ temp + pd, start = list(D = 1, g0 = 0.1, sigma = 9000), details = list(userdist = networkdistance, relativeD = TRUE), trace = FALSE, steptol = 1e-4) su.m2 <- secr.fit(capthist = GR.su, mask = list(mask2019, mask2020, mask2021), model = D ~ temp + temp2 + pd, 225 start = list(D = 1, g0 = 0.1, sigma = 9000), details = list(userdist = networkdistance, relativeD = TRUE), trace = FALSE, steptol = 1e-3) AIC selection of the best model kable(AIC(su.m0, su.m1, su.m2) %>% arrange(AICc)) model detectfn npar logLik AIC AICc dAICc AICcwt su.m2 D~temp + temp2 + pd g0~1 sigma~1 halfnormal 5 -1597.992 3205.983 3207.000 0.000 0.9743 su.m1 D~temp + pd g0~1 sigma~1 halfnormal 4 -1602.801 226 3213.601 3214.268 7.268 0.0257 su.m0 D~pd g0~1 sigma~1 halfnormal 3 -1611.034 3228.067 3228.461 21.461 0.0000 Predict the temperature covariate from the top model Create a dataframe of new data over which to predict the model covs19 <- covariates(mask2019) covs20 <- covariates(mask2020) covs21 <- covariates(mask2021) pd <- c(covs19$pd, covs20$pd, covs21$pd) newdat <- data.frame(temp = seq(min(covs19$temp), # coldest year max(covs21$temp), # warmest year length.out = 100), pd = mean(pd)) %>% # Set pd to a fixed value for prediction mutate(temp2 = temp^2) Predict the temperature covariate from the best model su.m2: cov.predict.temp <- predict(su.m2, newdata = newdat, type = "response", se.fit = TRUE, realname s = "D") preds.temp <- NULL for(i in 1:100){ 227 preds.temp <- rbind(preds.temp, cov.predict.temp[[i]]) } As the data were standardized for analyis, here we provide the non-transformed temperature variable from covs so that the x-axis of the plot is sensible. preds.temp$backtransformed <- seq(min(covs$temps19), # coldest year max(covs$temps21), # warmest year length.out = 100) Plot the response in Arctic grayling relative tag density to temperature: ggplot(preds.temp, aes(x = backtransformed, y = estimate)) + geom_line(linewidth = .7) + geom_line(aes(x = backtransformed, y = (estimate - SE.estimate)), linetype = 2, linewidth = .5) + geom_line(aes(x = backtransformed, y = (estimate + SE.estimate)), linetype = 2, linewidth = .5) + geom_ribbon(aes(x = backtransformed, ymin = (estimate - SE.estimate) , ymax = (estimate + SE.estimate)), fill = "gray50", alpha = 0.2) + labs(x = "Summer Temperature (°C)", y = "Activity centers per Rkm") + theme_bw() 228 Plotting model spatial outputs This section plots model predictions of the relative densities and activity centres of tagged Arctic grayling during each summer season. Data prep Get session-specific model outputs: m2.19 <- update(su.m2, capthist = GR.su[[1]], mask = mask2019) m2.20 <- update(su.m2, capthist = GR.su[[2]], mask = mask2020) m2.21 <- update(su.m2, capthist = GR.su[[3]], mask = mask2021) Here we represent activity centres for each tag as the pixel with the maximum probability of detection. These can be thought of as point estimates of average tag locations over the modeling window or the centre of home range selection. centers19 <- fxi.secr(m2.19) tags19 <- names(centers19) act.cent19 <- data.frame(tag = 1:length(centers19), x = NA, y = NA) for (i in 1:length(centers19)) { tmp2 <- max(centers19[[i]]) # Find maximum value in pdf for tag i tmp3 <- match(tmp2, centers19[[i]]) # Match value to index tmp4 <- mask2019[tmp3,] # Assign x-y coordinates of activity centre act.cent19[i,]$tag <- tags19[i] # Populate dataframe rows act.cent19[i,]$x <- tmp4$x act.cent19[i,]$y <- tmp4$y } centers20 <- fxi.secr(m2.20) tags20 <- names(centers20) act.cent20 <- data.frame(tag = 1:length(centers20), x = NA, y = NA) for (i in 1:length(centers20)) { tmp2 <- max(centers20[[i]]) # Find maximum value in pdf for tag i tmp3 <- match(tmp2, centers20[[i]]) # Match value to index tmp4 <- mask2020[tmp3,] # Assign x-y coordinates of activity centre act.cent20[i,]$tag <- tags20[i] # Populate dataframe rows act.cent20[i,]$x <- tmp4$x act.cent20[i,]$y <- tmp4$y } centers21 <- fxi.secr(m2.21) tags21 <- names(centers21) 229 act.cent21 <- data.frame(tag = 1:length(centers21), x = NA, y = NA) for (i in 1:length(centers21)) { tmp2 <- max(centers21[[i]]) # Find maximum value in pdf for tag i tmp3 <- match(tmp2, centers21[[i]]) # Match value to index tmp4 <- mask2021[tmp3,] # Assign x-y coordinates of activity centre act.cent21[i,]$tag <- tags21[i] # Populate dataframe rows act.cent21[i,]$x <- tmp4$x act.cent21[i,]$y <- tmp4$y } Create plots Working with relative densities in acoustic telemetry with multiple sessions may require some fine tuning of the color palette relative to the densities predicted in each session. The built-in plotting functions are useful for technical work and the spatial data can be saved as shapefiles for making further customized visualizations in a preferred GIS using the package sf (Pebesma and Bivand 2023, Pebesma 2018). Model outputs here differ slightly from what is presented in the manuscript as the bull trout covariate was not used in the models. cols19 <- terrain.colors(n = 13, rev = TRUE) plot(predictDsurface(m2.19), col = cols19, title = "Relative D") points(act.cent19$x, act.cent19$y) 230 cols20 <- terrain.colors(n = 7, rev = TRUE) plot(predictDsurface(m2.20), col = cols20, title = "Relative D") points(act.cent20$x, act.cent20$y) cols21 <- terrain.colors(n = 11, rev = TRUE) plot(predictDsurface(m2.21), col = cols21, title = "Relative D") points(act.cent21$x, act.cent21$y) 231 Save outputs for visualization in GIS Ensure that the GIS project is closed (or shapefiles are unloaded in the GIS environment) when running this code so that files can update. tag.cen19 <- act.cent19%>% st_as_sf(coords = c('x', 'y'), crs = 32610) rd.19 <- predictDsurface(m2.19) tag.den19 <- rd.19 %>% st_as_sf(coords = c('x', 'y'), crs = 32610) %>% mutate(D.0 = covariates(rd.19)$D.0) write_sf(tag.den19, 'outputs/GRtag.den2019.shp') write_sf(tag.cen19, 'outputs/GRtag.cen2019.shp') tag.cen20 <- act.cent20%>% st_as_sf(coords = c('x', 'y'), crs = 32610) rd.20 <- predictDsurface(m2.20) tag.den20 <- rd.20 %>% st_as_sf(coords = c('x', 'y'), crs = 32610) %>% mutate(D.0 = covariates(rd.20)$D.0) write_sf(tag.den20, 'outputs/GRtag.den2020.shp') write_sf(tag.cen20, 'outputs/GRtag.cen2020.shp') 232 tag.cen21 <- act.cent21%>% st_as_sf(coords = c('x', 'y'), crs = 32610) rd.21 <- predictDsurface(m2.21) tag.den21 <- rd.21 %>% st_as_sf(coords = c('x', 'y'), crs = 32610) %>% mutate(D.0 = covariates(rd.21)$D.0) write_sf(tag.den21, 'outputs/GRtag.den2021.shp') write_sf(tag.cen21, 'outputs/GRtag.cen2021.shp') 233 Appendix B. Expanded radio tag inspection plots 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261