Canadian Science Advisory Secretariat (CSAS) Research Document 2022/019 Arctic Region and Ontario and Prairie Region March 2022 Estimated population abundance and probability of population decline for Northern Hudson Bay narwhal (Monodon monoceros) Brooke A. Biddlecombe and Cortney A. Watt Arctic and Aquatic Research Division Fisheries and Oceans Canada 501 University Crescent Winnipeg, Manitoba, R3T 2N6 Foreword This series documents the scientific basis for the evaluation of aquatic resources and ecosystems in Canada. As such, it addresses the issues of the day in the time frames required and the documents it contains are not intended as definitive statements on the subjects addressed but rather as progress reports on ongoing investigations. Published by: Fisheries and Oceans Canada Canadian Science Advisory Secretariat 200 Kent Street Ottawa ON K1A 0E6 http://www.dfo-mpo.gc.ca/csas-sccs/ csas-sccs@dfo-mpo.gc.ca © Her Majesty the Queen in Right of Canada, 2022 ISSN 1919-5044 ISBN 978-0-660-42534-4 Cat. No. Fs70-5/2022-019E-PDF Correct citation for this publication: Biddlecombe, B.A., and Watt, C.A. 2022. Estimated population abundance and probability of population decline for Northern Hudson Bay narwhal (Monodon monoceros). DFO Can. Sci. Advis. Sec. Res. Doc. 2022/019. iv + 19 p. Aussi disponible en français : Biddlecombe, B.A., et Watt, C.A. 2022. Abondance estimée de la population et probabilité de déclin de la population pour le narval (Monodon monoceros) du nord de la baie d’Hudson. Secr. can. des avis sci. du MPO. Doc. de rech. 2022/019. iv + 20 p. http://www.dfo-mpo.gc.ca/csas-sccs/ mailto:csas-sccs@dfo-mpo.gc.ca iii TABLE OF CONTENTS ABSTRACT .................................................................................................................................. iv INTRODUCTION .......................................................................................................................... 1 METHODS .................................................................................................................................... 2 MODEL SPECIFICATION ......................................................................................................... 2 Priors ..................................................................................................................................... 2 Parameter estimation and model diagnostics ....................................................................... 3 Future projections and harvest scenarios ............................................................................. 4 RESULTS ..................................................................................................................................... 4 DISCUSSION ................................................................................................................................ 5 ACKNOWLEDGEMENTS ............................................................................................................. 6 REFERENCES CITED .................................................................................................................. 6 TABLES AND FIGURES ............................................................................................................... 9 APPENDIX .................................................................................................................................. 17 iv ABSTRACT Narwhal are an important subsistence harvest species for Inuit communities. The Northern Hudson Bay (NHB) population is spatially and genetically distinct from other populations of narwhal in Canada and Greenland. The NHB narwhal population has been assessed through periodic aerial surveys from 1981 – 2018. Survey estimates in 1982 and 2000 were negatively biased, as they did not factor in perception bias and were conducted and analyzed using different protocols. The 2011 and 2018 surveys addressed these biases. As a result, the 1982 and 2000 survey estimates were adjusted prior to population trend analysis using ratios that were calculated by re-analysing the 2011 survey data using methods similar to those applied in 1982 and 2000. To estimate the population trajectory and predict future population trends under various harvest scenarios, a Bayesian population model was fit to these four fully-adjusted aerial survey estimates and harvest data from 1951 – 2018. The preferred model resulted in a 2019 population estimate of 14,377 (95% CI 10,265 – 20,370), and an estimated starting population of 7,164 (95% CI 1,447 – 19,155) in 1951. Using a risk-based model approach to estimate the probability of population decline over 10 years, the model predicted a 0%, 20%, 40%, 50%, 60%, 80%, and 100% probability of decline with annual harvest quotas of 0, 63, 83, 93, 108, 173, and 450 narwhal per year. The incorporation of the updated 2018 survey estimate and the adjustment of past survey estimates provides greater confidence in the model results and, thus, in the risk-based approach for assessing harvest levels. The 2019 modelled estimated abundance was robust to input parameters in a series of model runs. Potential Biological Removal (PBR) was also calculated for the 2019 abundance estimate from the model and resulted in a removal threshold of 188 whales annually and a Total Allowable Landed Catch (TALC) of 151 narwhal. The risk-based approach considers the probability of population decline at different harvest levels, while the PBR inherently aims to keep the population at the maximum net productivity level (MNPL), which may result in a population decline if the population is already at or above MNPL. Management goals should be well defined to determine whether a PBR or a risk-based model approach is more appropriate for the NHB narwhal population. 1 INTRODUCTION Narwhal (Monodon monoceros) are distributed across both Canadian and Greenland Arctic seas. There are three recognized narwhal populations, Northern Hudson Bay (NHB), Baffin Bay (BB), and East Greenland (EG). NHB narwhal summer in northern Hudson Bay and migrate through Hudson Strait to spend winters in southern Davis Strait (Richard 1991, Westdal et al. 2010) (Figure 1). NHB narwhal are spatially and genetically distinct from other narwhal populations in Canada and Greenland (Petersen et al. 2011). Narwhals are a culturally important species for Inuit communities, with a long history of subsistence harvest (Richard and Pike 1993). Non-commercial harvests had sparse records until the creation of a harvest quota in 1977 (Stewart 2008). Subsistence harvest of NHB narwhal occurs primarily near the community of Naujaat (formerly Repulse Bay), and to lesser extents in Salluit, Kinngait, Igloolik, Chesterfield Inlet, Coral Harbour, Kimmirut, Rankin Inlet, and Whale Cove. NHB narwhal have a minimal commercial harvest history, as bowhead whales were the focus of commercial whaling in this region (Stewart 2008). To ensure a sustainable subsistence harvest, surveys and resulting abundance estimates need to be consistently updated. Aerial surveys during the summer months are used to estimate narwhal abundance, as narwhal typically show long-term site fidelity to summering regions (Richard 1991, Heide-Jørgensen et al. 2003, Westdal et al. 2010). Surveys of the NHB population were flown in each year from 1981 – 1984, and again in 2000, 2011, and 2018 (Richard 1991, Bourassa 2003, Asselin et al. 2012, Watt et al. 2020). The methodology of the 1980’s surveys varied between reconnaissance, random, or systematic in design, and using visual observers or photographic count methods. All later surveys used systematic visual observation methods. Thus, of the early surveys, only the estimate from 1982 is appropriate for use in population modelling because it was the only systematic visual survey during that time (Asselin and Ferguson 2013). Despite being systematic in design, the 1982 survey did not account for perception bias and was analyzed as a strip transect, which assumed that all narwhals within a 600 m strip width were observed. However, detectability actually declines with increasing distance from the aircraft (Buckland et al. 2001) and, therefore, this estimate is negatively biased. The survey conducted in 2000 did not account for perception bias, also resulting in negative bias. To make the 2011 survey, which was analyzed using mark-recapture distance analysis and, therefore, incorporated perception bias, comparable with the 1982 and 2000 survey, the 2011 survey was reanalyzed using the methods of the two older surveys. In this way a ratio was calculated by which the 1982 and 2000 surveys could be adjusted to approximately equate them to the methods used in 2011. The ratio used to adjust the 1982 survey was 2.56, and the 2000 survey ratio was 2.29 (Asselin and Ferguson 2013). The NHB narwhal population abundance has been modelled previously by Kingsley et al. (2013), who adjusted the 2011 survey estimate. However, Kingsley et al.’s (2013) model estimates had a high degree of associated uncertainty, which made a risk-based approach for assessing harvest levels unreliable until future surveys were performed (Kingsley et al. 2013). As such, a new aerial survey was completed in the summer of 2018 (Watt et al. 2020), which allows us to update the population model, potentially providing increased certainty in the ability of the model to assess the impact of various harvest levels. The objective of this study is to use harvest data and survey estimates to model the population trends of NHB narwhal and predict future trajectories under various harvest scenarios. 2 METHODS MODEL SPECIFICATION Bayesian Markov Chain Monte Carlo (MCMC) methods were used to fit a stochastic stock production model with density dependence acting on the population growth rate. This model utilizes aerial survey and harvest data to estimate population dynamics. To separate two stochastic processes, observation process and state process, a hierarchical state-space model was developed. Observation process describes observation error, which results from data collection and estimation of abundance. State process describes process error, which reflects natural variability in population dynamics (de Valpine and Hastings 2002). The state process was defined using a discrete formulation of the Pella and Tomlinson model (Pella and Tomlinson 1969, Innes and Stewart 2002), which describes the change in true population size over time. The population size at time t is a multiple of the previous year’s population with harvest removals deducted. 𝑁𝑁𝑡𝑡 = 𝑁𝑁𝑡𝑡−1 ∗ �1 + (R𝑚𝑚𝑚𝑚𝑚𝑚) ∗ �1 − � 𝑁𝑁𝑡𝑡−1 𝐾𝐾 � θ �� ∗ ε𝑝𝑝 − 𝑅𝑅𝑡𝑡 Where: Rmax is the maximum growth rate or rate of population increase, 𝐾𝐾 is the environmental carrying capacity, 𝜃𝜃 defines the shape of the density-dependent function (theta), 𝜀𝜀𝑝𝑝 is a stochastic term for the process error, 𝑅𝑅t are the harvest removals for that year, calculated as reported catches, 𝐶𝐶t , that are corrected for the proportion of animals that were struck and lost, 𝑆𝑆&𝐿𝐿: 𝑅𝑅𝑡𝑡 = 𝐶𝐶𝑡𝑡 ∙ (1 + 𝑆𝑆&𝐿𝐿) The observation process links observed data to true population size. We determined survey error by using a multiplicative error (𝜀𝜀𝑆𝑆𝑡𝑡) term to link true population size (𝑁𝑁𝑡𝑡) to aerial survey estimates (𝑆𝑆𝑡𝑡). ln( 𝑆𝑆𝑡𝑡) = ln(𝑁𝑁𝑡𝑡) + 𝜀𝜀𝑆𝑆𝑡𝑡 The model was fit using harvest data from 1951 – 2018 (Table 1) and aerial survey data from 1982, 2000, 2011, and 2018 (Table 2) to estimate population dynamics for 1951 – 2018. Harvest data was available starting in 1911, but these data were too temporally sparse before 1951 to be included in the model (Stewart 2008). Surveys from 1982 and 2000 were adjusted to account for differences in survey design and analysis, and to incorporate an adjustment for perception bias, by multiplying them by 2.29 and 2.56, respectively (Asselin and Ferguson 2013; Table 2). The 2011 and 2018 survey estimates had been previously adjusted for perception bias. All survey estimates were adjusted to account for animals that were diving and not visible during the survey (availability bias). Priors Prior distributions for the random variables included in the model were informed by traditional knowledge, previous cetacean population models (i.e., Wade 1998, Marcoux and Hammill 2016), and initial exploratory runs of our model. Maximum population growth rate (Rmax) for cetaceans is often assumed to be 0.04 (Wade 1998), therefore we fixed Rmax at 0.04 for some initial model runs but allowed the final model to estimate Rmax using a uniform distribution ranging from 0.01 – 0.07 (Table 3). The shaping parameter θ (theta) was explored with multiple 3 priors in initial runs, from which we determined our final model run was allowed to vary between one and five (Table 3). Starting population size was assigned a uniform prior distribution from 1,000 to 20,000 (Table 3). Carrying capacity (K) was assigned a uniform distribution with a minimum of 10,000 and a maximum of 30,000 (Table 3). Harvest data underestimate the true number of narwhal removals because some animals are wounded or killed but cannot be recovered by hunters and are therefore referred to as struck and lost (S&L). Naujaat is the primary harvest community for NHB narwhal, thus, we based our S&L prior on the loss rate correction values for Naujaat in 1999 – 2005 (Richard 2008). Mean loss rate correction in Naujaat was 0.25, which we used to assign a moderately informative prior following a beta(3.5, 10) distribution, with a median of 0.25, and quartiles of 0.17 and 0.33. The stochastic process error term had a log-normal distribution with a zero location parameter. The precision parameter for process error was assigned an informative gamma(1.5, 0.00001) distribution. The resulting process error term had quartiles of 0.998 and 1.002. The narrow stochastic process range reflects the assumption that because narwhal are a long-lived species, their stock dynamics have low inter-annual variability. Uncertainty associated with aerial surveys was incorporated into model fitting by informing the prior distribution for survey error. The formulation of the prior for the survey error was weighted by the standard error (SE) for each of the abundance estimates from the four surveys incorporated in the model. The survey error term was given a lognormal distribution with a zero location parameter. The precision parameter for survey error was specific to each survey year using a precision value calculated from each SE independently. Parameter estimation and model diagnostics A Gibbs sampler algorithm implemented in JAGS (Plummer 2003) was used to obtain posterior estimates for the parameters included in the model. Results were examined in RStudio (R Core Team 2020, version 1.2.5033) using packages R2jags (Su and Yajima 2020, version 0.6-1) and coda (Plummer et al. 2006, version 0.19-3). We used a MCMC chain simulation method which required us to check the convergence of sample values in each parameter to a stationary distribution. Model convergence, mixing, and autocorrelation within parameters were investigated in initial runs of the model code. Geweke’s Diagnostic was used to test similarity between different sections of the chains and Geweke plots were visually inspected to test for mixing within each chain (Geweke 1996). Convergence between chains was validated by comparing the width of the 80% credible interval of the chains pooled with the mean widths of the 80% credible interval for each of the chains individually using the Brooks-Gelman-Rubin (BGR) diagnostic test (Brooks and Gelman 1998). Sensitivity of model results to our updated prior for survey error precision informed by the SE from surveys was evaluated by comparing resulting population trajectories and whether the widths of confidence intervals around the median reflected the SE from each survey. This was compared to initial runs using a gamma(2.5, 0.4) prior distribution which is used in previous population models for Arctic toothed whales (Marcoux and Hammill 2016). Model variations were also run allowing θ to be estimated within the uniform prior described above and with θ = 1, to compare the effects of linear and non-linear density dependence. All initial model runs used three chains and had the number of iterations corrected to equal 10,000 after a burn in of 6,000 and thinning value of 30. After assessing outputs from all initial runs, we determined a preferred model which was the most parsimonious and biologically relevant model (Table 3). The preferred model was run using five chains, and 100,000 iterations after a burn in of 60,000, with thinning maintained at 30. 4 Future projections and harvest scenarios To predict stock trajectory and sustainable yield of this population, the model was extended 10 years into the future under 10 different harvest scenarios. A relatively short period of 10 years was used due to the frequency with which harvest quotas are re-assessed for marine mammals. Recent survey estimates suggest that the population is increasing, so our harvest scenarios ranged from 0 – 450 narwhals taken annually, in intervals of 50. The impact of each harvest level was assessed by calculating the probability of population decline over 10 years. This is the risk-based approach for assessing harvest levels and is used for populations that are considered data rich. To further explore results, we estimated Potential Biological Removal (PBR). Unlike the risk- based approach, PBR does not focus on probabilities of increase or decline, but rather its built- in management objective is to quantify the maximum number of animals that can be removed from a marine mammal stock while still allowing the stock to reach or stay at its Maximum Net Productivity Level (MNPL) after 100 years (Wade 1998). PBR is the default method for estimating sustainable yield for stocks that are considered data poor. PBR was calculated using a recovery factor (FR) of 0.75, since this population is abundant, but data limited (Hammill et al. 2017). The equation for PBR threshold is: 𝑃𝑃𝑃𝑃𝑅𝑅 = 𝑁𝑁𝑚𝑚𝑚𝑚𝑚𝑚 ∙ 0.5 ∙ 𝑅𝑅𝑚𝑚𝑚𝑚𝑚𝑚 ∙ 𝐹𝐹𝑅𝑅 Where: 𝑅𝑅max is the maximum rate of population increase. The default value for cetaceans is 0.04 𝐹𝐹R is a recovery factor (between 0.1 and 1), and 𝑁𝑁min is the estimated population size using the 20th-percentile of the posterior distribution resulting from the model or the 20th-percentile of the log-normal distribution (Wade 1998) of the aerial survey estimate. As explained below, we used the 2019 population estimate from the model to calculate the Nmin value. PBR estimates total removals from a population, which includes animals that are harvested, those that are caught but not landed, unreported harvests, and other sources of human caused mortality. Therefore, the Total Allowable Landed Catch (TALC) is: 𝑇𝑇𝑇𝑇𝐿𝐿𝐶𝐶 = 𝑃𝑃𝑃𝑃𝑅𝑅/𝐿𝐿𝑅𝑅𝐶𝐶 where LRC is the loss rate correction. LRC is estimated from the S&L value from model results to calculate TALC with modelled 2019 abundance estimates. RESULTS Variables carrying capacity (𝐾𝐾), population size in 2019 (N2019), process error, initial population size (N1951), and S&L rate showed rapid convergence for all chains. Autocorrelation was low in all variables, and the BGR statistics were close to 1 for all variables, indicating convergence of each chain and each model (see Appendix). In early model considerations, multiple runs with varied priors were performed and the model was robust to changes in priors. The preferred model is presented here. The model shows a gradual increase in population abundance over time, which is consistent with the aerial survey estimates. The model posterior distributions resulted in a median maximum growth rate (Rmax) of 0.038, the median carrying capacity (K) was estimated as 16,779, and θ was estimated as 2.94. The median starting population in 1951 was estimated as 5 7,164 (95% CI 1,447 – 19,155), and N2019 was estimated as 14,377 (95% CI 10,265 – 20,370) (Figure 2). For a risk-based approach to assessing harvest levels, the population trajectories for 10 harvest scenarios were plotted using the model results (5 of 10 harvest levels shown in Figure 3). These trajectories determined that after 10 years, the probability of population decline was 10% for an annual catch level of 54, 30% for a catch level of 73, 50% for a catch level of 93, 70% for a catch level of 173, and 90% for a catch level of 243 narwhal (Table 4, Figure 4). The 2019 modelled population estimate from our most parsimonious and biologically relevant model was used in the PBR calculation because there was more confidence in the model estimate than in the 2018 aerial survey estimate. Based on the 2019 population abundance estimate, a recovery factor FR of 0.75, and a default maximum growth rate of 4% per year, the PBR threshold was 188 narwhals. TALC was calculated as 151, using the model estimate for S&L (0.245) and the 2019 abundance estimate. DISCUSSION In this study, we used Bayesian methodology to fit a stochastic Pella-Tomlinson model (Pella and Tomlinson 1969, Innes and Stewart 2002) to narwhal aerial survey data (1982 – 2018) and a time series of catch history data (1951 – 2018). The resulting model had a shaping parameter (θ) greater than one, reflecting a convex growth response (i.e., density feedback occurs at higher proportions of K). The estimated median population started at 7,164 in 1951 and increased over time to a 2019 estimate of 14,377 narwhals. Narwhal maximum population growth can vary, depending on the vital rates of the population (Kingsley 1989). In our model Rmax was estimated as 0.038 which agrees with the commonly assumed Rmax of 0.04 for cetaceans (lambdamax = 1.04) (Wade 1998). Adjustments to the aerial survey results from 1982 and 2000 were necessary for consistency across the time series. Using adjusted survey estimates allowed the model to re-create a plausible population trajectory which shows a gradual increase over the time series. Allowing the model to estimate density dependence resulted in the population trajectory fitting the aerial survey estimates slightly better than when θ = 1. Marine mammals commonly exhibit non-linear density dependence, resulting in shaping parameters that fall between 1 and 7 (Taylor and DeMaster 1993) so the θ estimate of 2.94 is biologically plausible for narwhal as a long-lived species. The θ estimate translates to the population reaching maximum productivity in the range of 50 – 70% of K with growth slowing down after that point. The 2019 population estimate of 14,377 is 86% of the estimated K of 16,779, which suggests that the population has reached a point of reduced growth, as is evident by the flattened curve in the latter years of the time series. Aerial surveys from 1982 and 2000 were adjusted post hoc, which is a potential source of uncertainty in this model. Adjustments for these earlier surveys reflect both perception bias and changes in the analytical methods used for later surveys (Asselin and Ferguson 2013). However, there is potential that the adjustments result in over- or under-estimates of the true population abundance. A further artifact of using past aerial survey abundance estimates is that the latter two surveys (2011 and 2018) covered additional areas (namely Roes Welcome Sound and Wager Bay, which had surface abundance estimates of ~ 1100 and ~ 500 in 2011 and 2018, respectively) than the 1982 and 2000 surveys (Asselin et al. 2012, Watt et al. 2020). The two earlier surveys may have a negative bias due to their lack of coverage in regions where narwhals were present in 2011 and 2018. We could not correct for this potential bias as the timeline for any distribution changes in NHB narwhal is unknown. The types of uncertainties described are not uncommon for models that utilize historical data. 6 The model shows that the NHB narwhal population has increased over the entire time series from 1951 – 2018. There was a moderate increase in harvest levels in the early 2000’s which appears to have slightly affected the population trajectory, although overall growth still occurred. Using PBR, the TALC value was calculated as 151. Using the risk-based approach, a harvest of 151 reflects a ~ 76% probability of population decline in 10 years. A 5% probability of decline under the risk-based approach is achieved at a harvest level of 36. The risk-based method assesses the probability of population decline from the current abundance whereas PBR strives to keep a population at or above its MNPL, which can allow for decline if the population is already above that level. As compared to previous NHB narwhal population models (Kingsley et al. 2013), there is more certainty associated with the model results from this study as a result of the addition of a more recent survey estimate and the incorporation of adjustments to earlier survey estimates. As a consequence of the updates applied to this model, there is more confidence in the current population estimate which can be used to assess harvest levels either using population modelling with the risk-based approach or via PBR. The preferred approach for informing the NHB narwhal harvest quota will depend on management objectives. ACKNOWLEDGEMENTS We thank the Naujaat Hunters and Trappers Association and J. Young for providing harvest information. We thank the late M. Kingsley for developing the initial population model, T. Doniol- Valcroze, A. Mosnier, and R. Hobbs for model improvements, and M. Hammill for modelling advice. REFERENCES CITED Asselin, N.C., Ferguson, S.H., Richard, P.R. and Barber, D.G. 2012. Results of narwhal (Monodon monoceros) aerial surveys in northern Hudson Bay, August 2011. DFO Can. Sci. Advis. Sec. Res. Doc. 2012/037. iii + 23 p. Asselin, N.C., and Ferguson, S.H. 2013. A re-analysis of northern Hudson Bay narwhal surveys conducted in 1982, 2000, and 2011. DFO Can. Sci. Advis. Sec. Res. Doc. 2013/019. iv + 9 p. Bourassa, M.N. 2003. Inventaires de la population de narvals (Monodon monoceros) du nord de la baie d’Hudson et analyse des changements démographiques depuis 1983. Thesis (M.Sc) Université du Québec à Rimouski, Rimouski, QC. xii + 69 p. Brooks, S.P., and Gelman, A. 1998. Alternative methods for monitoring convergence of iterative simulations. J. Comput. Graph. Stat. 7: 434–455. Buckland, S.T., Anderson, D.R., Burnham, K.P., Laake, J.L., Borchers, D.L., and Thomas, L. 2001. Introduction to distance sampling: estimating abundance of biological populations. Oxford University Press, Oxford, UK. 448 p. de Valpine, P., and Hastings, A. 2002. Fitting population models incorporating process noise and observation error. Ecol. Monogr. 72(1): 57–76. doi:10.2307/3100085. Geweke, J. 1996. Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. In Bayesian Statistics 4. Edited by J.M. Bernardo, J.M. Berger, A.P. Dawid and A.F.M.. Oxford University Press, Oxford, UK. pp. 169–193. Hammill, M.O., Stenson, G.B., and Doniol-Valcroze, T. 2017. A management framework for Nunavik beluga. DFO Can. Sci. Advis. Sec. Res. Doc. 2017/060. v + 34 p. https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2012/2012_037-eng.html https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2012/2012_037-eng.html https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2013/2013_019-eng.html https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2013/2013_019-eng.html https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2017/2017_060-eng.html https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2017/2017_060-eng.html 7 Heide-Jørgensen, M.P., Dietz, R., Laidre, K.L., Richard, P., Orr, J., and Schmidt, H.C. 2003. The migratory behaviour of narwhals (Monodon monoceros). Can. J. Zool. 81(8): 1298–1305. Innes, S., and Stewart, R.E.A. 2002. Population size and yield of Baffin Bay beluga (Delphinapterus leucas) stocks. NAMMCO Sci. Publ. 4: 225–238. doi:10.7557/3.2846. Kemper, J.B. 1980. History of use of narwhal and beluga by Inuit in the Canadian eastern Arctic including changes in hunting methods and regulations. Rep. Intern. Whal. Comm. 30: 481–492. Kingsley, M. 1989. Population dynamics of the narwhal Monodon monoceros: an initial assessment (Odontoceti: Monodontidae). J. Zool. 219(2): 201–208. Kingsley, M.C.S., Asselin, N.C., and Ferguson, S.H. 2013. Updated stock-dynamic model for the Northern Hudson Bay narwhal population based on 1982-2011 aerial surveys. DFO Can. Sci. Advis. Sec. Res. Doc. 2013/011. v + 19 p. Mansfield, A.W., Smith, T.G., and Beck, B. 1975. The narwhal, Monodon monoceros, in eastern Canadian waters. J. Fish. Res. Board Can. 32(7): 1041–1046. Marcoux, M., and Hammill, M.O. 2016. Model estimates of Cumberland Sound beluga (Delphinapterus leucas) population size and total allowable removals. DFO Can. Sci. Advis. Sec. Res. Doc. 2016/077. iv + 35 p. Pella, J.J., and Tomlinson, P.K. 1969. A generalized stock production model. Inter-Am. Trop. Tuna Comm. 13(3): 420–496. Petersen, S.D., Tenkula, D., and Ferguson, S.H. 2011. Population genetic structure of narwhal (Monodon monoceros). DFO Can. Sci. Advis. Sec. Res. Doc. 2011/021. vi + 20 p. Plummer, M. 2003. A program for analysis of Bayesian graphical models using Gibbs sampling. In Proceedings of the 3rd International Workshop on Distributed Statistical Computing. Edited by K. Hornik, F. Leisch, and A. Zeileis. Technische Universität Wien, Vienna, Austria. pp. 1–10. Plummer, M., Best, N., Cowles, K., and Vines, K. 2006. CODA: Convergence Diagnosis and Output Analysis for MCMC. R News 6: 7–11. R Core Team. 2020. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Richard, P.R. 1991. Abundance and distribution of narwhals (Monodon monoceros) in northern Hudson Bay. Can. J. Fish. Aquat. Sci. 48: 276–283. Richard, P.R. 2008. On determining the Total Allowable Catch for Nunavut odontocete stocks. DFO Can. Sci. Advis. Sec. Res. Doc. 2008/022. iv + 12 p. Richard, P.R. and Pike, D.G. 1993. Small whale co-management in the eastern Canadian Arctic: a case history and analysis. Arctic 46(2): 138–143. Spiegelhalter, D.J., Best, N.G., Carlin, B.P., and van der Linde, A. 2002. Bayesian measures of model complexity and fit. J. Royal Statist. Soc. Ser. B Statist. Methodol. 64(4): 583–639. doi:10.1111/1467-9868.00353. Stewart, D.B. 2008. Commercial and subsistence harvests of narwhals (Monodon monoceros) from the Canadian eastern Arctic. Arctic Biological Consultants, Winnipeg, MB for Canada Department of Fisheries and Oceans, Winnipeg, MB. ii + 97 p. https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2013/2013_011-eng.html https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2013/2013_011-eng.html https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2016/2016_077-eng.html https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2016/2016_077-eng.html https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2011/2011_021-eng.html https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2011/2011_021-eng.html https://www.r-project.org/ https://www.dfo-mpo.gc.ca/csas-sccs/publications/resdocs-docrech/2008/2008_022-eng.htm 8 Strong, J.T. 1989. Reported harvests of narwhal, beluga and walrus in the Northwest Territories, 1948-1987. Can. Data Rep. Fish. Aquat. Sci. 734: iv + 14 p. Su, Y., and Yajima, M. 2020. R2jags: Using R to Run 'JAGS'. R package version 0.6-1. Taylor, B.L., and DeMaster, D.P. 1993. Implications of non-linear density dependence. Mar. Mammal Sci. 9(4): 360–371. doi:10.1111/j.1748-7692.1993.tb00469.x. Wade, P.R. 1998. Calculating limits to the allowable human-caused mortality of cetaceans and pinnipeds. Mar. Mammal Sci. 14(1): 1–37. Watt, C.A., Hornby, C., and Hudson, J. 2020. Narwhal (Monodon monoceros) abundance estimate from the 2018 aerial survey of the Northern Hudson Bay population. DFO Can. Sci. Advis. Sec. Res. Doc. 2020/073. iv + 15 p. Westdal, K. H., P. R. Richard, and J. R. Orr. 2010. Migration route and seasonal home range of the northern Hudson Bay narwhal (Monodon monoceros). In A Little Less Arctic: Top predators in the world’s largest northern inland sea, Hudson Bay. Edited by S.H. Ferguson, L.L Loseto, and M.L. Mallory. Springer, Dordrecht, Netherlands. pp. 71–92. https://publications.gc.ca/site/eng/9.562737/publication.html https://publications.gc.ca/site/eng/9.562737/publication.html https://cran.r-project.org/package=R2jags https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2020/2020_073-eng.html https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2020/2020_073-eng.html 9 TABLES AND FIGURES Table 1. Reported harvests for Northern Hudson Bay narwhal from 1951 – 2018. Data for 1951 – 1953 are from Stewart (2008), data for 1953 – 2018 are from Mansfield (1975), Kemper (1980), Strong (1989), and DFO harvest statistics 1988 – 2018 (unpublished data from DFO Resource Management. Year Reported catch Year Reported catch 1951 0 1985 27 1952 0 1986 16 1953 0 1987 7 1954 0 1988 16 1955 0 1989 38 1956 0 1990 16 1957 175 1991 17 1958 35 1992 36 1959 15 1993 20 1960 0 1994 14 1961 0 1995 6 1962 0 1996 14 1963 0 1997 45 1964 0 1998 0 1965 23 1999 28 1966 100 2000 157 1967 73 2001 45 1968 2 2002 108 1969 0 2003 67 1970 0 2004 43 1971 5 2005 120 1972 14 2006 86 1973 1 2007 94 1974 0 2008 88 1975 0 2009 29 1976 8 2010 118 1977 0 2011 103 1978 6 2012 92 1979 31 2013 56 1980 1 2014 115 1981 26 2015 96 1982 41 2016 47 1983 23 2017 73 1984 11 2018 113 10 Table 2. Estimates of Northern Hudson Bay narwhal abundance from aerial surveys in their summering range. Total Estimates are the values calculated from the surface estimates using an availability bias adjustment factor (Ca) of 2.8, and the coefficient of variation (CV) of the total estimate. Fully Adjusted Estimates have been adjusted for both perception and availability bias. CDS = conventional distance sampling, MRDS = mark-recapture distance sampling. Survey Year Total Estimate CV Observation Method Analysis Method Adjustment Ratio Fully Adjusted Estimate Citation 1981 n/a n/a Reconnaissance visual - - - Richard 1991 *1982 2,906 0.52 Systematic visual 600 m Strip 2.56 7,440 Richard 1991, Asselin and Ferguson 2013 1983 4,248 0.36 Systematic photographic - - - Richard 1991 1984 3,794 0.31 Systematic photographic - - - Richard 1991 *2000 4,978 0.40 Systematic visual CDS 2.29 11,401 Bourassa 2003 (CV in Richard 2008), Asselin and Ferguson 2013 *2011 12,485 0.26 Systematic visual MRDS 1.00 12,485 Asselin et al. 2012 *2018 19,232 0.28 Systematic visual MRDS 1.00 19,232 Watt et al. 2020 * denotes survey data that were used in the population dynamics model. 11 Table 3. Prior distributions, parameters, and hyper-parameters used in the population model. “dist.” denotes a hyper-parameter with its own prior distribution. Parameters Notation Prior distribution Hyper-parameters Values Survey error (t) εst Log-normal μs 0 τs dist. Precision (survey) τs * - - Process error (t) εpt Log-normal μp 0 τp dist. Precision (process) τp Gamma αp 1.5 βp 0.00001 Density dependence shape function θ Uniform Nupp 5 Nlow 1 Struck-and-lost S&L Beta αsl 3.5 βsl 10 Initial population N1951 Uniform Nupp 20,000 Nlow 1,000 Carrying capacity K Uniform Nupp 30,000 Nlow 10,000 Max. growth rate Rmax Uniform Nupp 0.07 Nlow 0.01 * weighted by the standard error from each survey 12 Table 4. Probability (P) that the NHB narwhal population, subjected to different levels of annual landed catch, will decline from the modelled 2019 population abundance estimate after 10 years of harvest. P (%) Landed Catch 0 0 10 54 20 63 30 73 40 83 50 93 60 108 70 135 80 173 90 243 100 450 13 Figure 1. Map indicating the three regions in Nunavut, and the summering area for the Northern Hudson Bay narwhal population (dark blue). 14 Figure 2. Model estimates of NHB narwhal abundance from the model fitted to adjusted aerial survey estimates from 1982, 2000, 2011 and 2018 (black dots ± 95% log-normal CI) and harvest data from 1951–2018. Solid line shows the median estimates and dashed lines show the 2.5th, 25th, 75th, and 97.5th quantiles. 15 Figure 3. Future population projections for NHB narwhal under 5 different harvest scenarios based on population estimates from the model, fitted to adjusted aerial survey estimates from 1982, 2000, 2011 and 2018 (black dots ± 95% log-normal CI). Solid line is the median abundance estimate and the dashed lines are the 95% confidence intervals. Vertical dotted line shows separation between the modelled abundance based on aerial survey estimates and the future projections. 16 Figure 4. Probability of the NHB narwhal stock decreasing from the 2019 abundance estimate after 10 years of harvest, estimated by the model, as a function of the number of reported narwhal removed from the stock each year. Dotted lines show the corresponding probability of decline (y-axis) for various annual levels of harvest (x-axis). 17 APPENDIX Detailed outputs from the model to examine population abundance trends of Northern Hudson Bay narwhal using catch history data from 1951 – 2018, and aerial survey estimates from 1982, 2000, 2011, and 2018. Table A1. Model outputs for NHB narwhal stock using 1951–2018 catch history and 1982–2018 adjusted survey estimates. The mean, standard deviation (SD), 2.5th , 25th, 50th, 75th and 97.5th quantiles are given for the following model parameters and their priors: carrying capacity (K), shaping parameter (theta), process error (prec.process), survey precision (prec.survey), starting population (startpop), struck and lost (S&L), and population size in 2019 (N2019). 𝑅𝑅� is the Brooks-Gelman-Rubin statistic; values near one indicate convergence of chains. N.eff is the number of effective runs after considering autocorrelation. Mean sd 2.50% 25% 50% 75% 97.50% Rhat n.eff K 18044.534 4846.625 11377.84 14324.234 16779.243 21070.311 28935.97 1.001 5.00E+05 K.prior 19994.333 5770.36 10499.53 14991.805 19994.798 24980.759 29499.297 1.001 5.00E+05 theta 2.957 1.16 1.092 1.943 2.937 3.961 4.895 1.001 5.00E+05 theta.prior 2.999 1.154 1.101 2 2.999 3.999 4.899 1.001 5.00E+05 deviance 77.473 1.765 74.764 76.162 77.52 78.302 81.604 1.001 5.00E+05 Rmax 0.039 0.016 0.012 0.025 0.038 0.052 0.068 1.001 5.00E+05 Rmax.prior 0.04 0.017 0.012 0.025 0.04 0.055 0.069 1.001 3.00E+05 prec.process 150011.35 122358.2 10759.791 60701.106 118449.909 205465.538 466990.681 1.001 38000 prec.process.prior 149919.295 122634.137 10621.624 60362.222 118148.845 205384.57 467770.558 1.001 5.00E+05 prec.surv 10.585 5.502 2.693 6.555 9.644 13.597 23.793 1.001 370000 prec.surv.prior 10.561 5.495 2.69 6.53 9.624 13.591 23.719 1.001 5.00E+05 startpop 8341.686 5241.921 1446.671 3882.189 7163.573 12199.53 19154.801 1.001 5.00E+05 startpop.prior 10510.619 5484.469 1479.906 5757.679 10514.792 15270.275 19525.115 1.001 5.00E+05 struck.and.lost 0.257 0.114 0.072 0.171 0.245 0.33 0.509 1.001 5.00E+05 struck.and.lost.prior 0.259 0.115 0.073 0.173 0.247 0.333 0.513 1.001 5.00E+05 N2019 14622.06 2592.588 10264.86 12757.33 14376.96 16222.03 20370.09 1.001 5.00E+05 18 Figure A1. The model was fit to harvest data (1951 – 2018) and adjusted aerial survey data (1982 – 2018). Plots show change in A) autocorrelation, B) cross correlation, and C) BGR test results. 19 Figure A2. Priors (lines) and posteriors (histograms) for A) carrying capacity, B) initial population, C) S&L, D) Rmax, and E) theta. ABSTRACT INTRODUCTION METHODS MODEL SPECIFICATION Priors Parameter estimation and model diagnostics Future projections and harvest scenarios RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES CITED TABLES AND FIGURES APPENDIX