Canadian Science Advisory Secretariat (CSAS) Research Document 2023/087 Ontario and Prairie Region December 2023 Recovery Potential Modelling of Purple Wartyback (Cyclonaias tuberculata) in Canada Adam S. van der Lee and Marten A. Koops Fisheries and Oceans Canada Great Lakes Laboratory for Fisheries and Aquatic Sciences 867 Lakeshore Rd. Burlington ON L7S 1A1 Canada 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 © His Majesty the King in Right of Canada, as represented by the Minister of the Department of Fisheries and Oceans, 2023 ISSN 1919-5044 ISBN 978-0-660-68881-7 Cat. No. Fs70-5/2023-087E-PDF Correct citation for this publication: van der Lee, A.S. and Koops, M.A. 2023. Recovery Potential Modelling of Purple Wartyback (Cyclonaias tuberculata) in Canada. DFO Can. Sci. Advis. Sec. Res. Doc. 2023/087. iv + 24 p. Aussi disponible en français : van der Lee, A.S. et Koops, M.A. 2023. Modélisation du potentiel de rétablissement de la mulette verruqueuse (Cyclonaias tuberculata) au Canada. Secr. can. des avis sci. du MPO. Doc. de rech. 2023/087. iv + 28 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 .................................................................................................................................... 1 SOURCES ................................................................................................................................ 1 THE POPULATION MODEL ..................................................................................................... 2 PARAMETERIZATION .............................................................................................................. 3 Life-History ............................................................................................................................ 4 Density-Dependence ............................................................................................................. 7 Stochasticity .......................................................................................................................... 8 IMPACT OF HARM ................................................................................................................... 8 Elasticity of λ ......................................................................................................................... 9 Simulation ............................................................................................................................. 9 RECOVERY TARGETS ............................................................................................................ 9 Abundance: Minimum Viable Population (MVP) ................................................................... 9 Habitat: Minimum Area for Population Viability (MAPV) ..................................................... 10 RECOVERY TIMES ................................................................................................................ 11 RESULTS ................................................................................................................................... 11 IMPACT OF HARM ................................................................................................................. 11 Elasticity of λ ....................................................................................................................... 11 Simulation ........................................................................................................................... 14 RECOVERY TARGETS .......................................................................................................... 14 Abundance: Minimum Viable Population (MVP) ................................................................. 14 Habitat: Minimum Area for Population Viability (MAPV) ..................................................... 16 RECOVERY TIMES ................................................................................................................ 16 DISCUSSION .............................................................................................................................. 17 UNCERTAINTIES ................................................................................................................... 19 ELEMENTS ............................................................................................................................. 21 REFERENCES CITED ................................................................................................................ 22 iv ABSTRACT The Committee on the Status of Endangered Wildlife in Canada (COSEWIC) assessed Purple Wartyback (PWB, Cyclonaias tuberculata) as Threatened in Canada. Here population modelling is presented to assess the impacts of harm, determine abundance and habitat recovery targets, and conduct long-term projections of population recovery in support of a recovery potential assessment (RPA). The model incorporated parameter uncertainty, environmental stochasticity, and density-dependence into population projections. The analysis demonstrated that PWB populations were most sensitive to perturbations to adult survival under most circumstances. As population growth rate (λ) increased the sensitivity of juvenile survival to perturbation increased and surpassed adult survival sensitivity when λ > 1.2. Estimates of the level of harm that would reduce population growth rate to 1 were estimated for populations in the Sydenham and Thames rivers. Population viability analysis was used to identify potential recovery targets. Demographic sustainability, (i.e., a self-sustaining population over 250 years) can be achieved with population sizes of ~2,800 (CI: 1,900–4,000) adults. It was estimated that populations of minimum viable population (MVP) size would require 623.3 m2 (CI: 251.9–1,396.9) and 2,900 m2 (CI: 301.5–17,166.3) in the Sydenham and Thames rivers respectively. Therefore, there is sufficient habitat to support PWB populations in both systems. 1 INTRODUCTION Purple Wartyback (PWB, Cyclonaias tuberculata) is a medium sized freshwater mussel (Unionidae) endemic to eastern North America. In Canada, PWB was historically found in Lake Erie, and the Detroit, Sydenham, Thames and Ausable rivers. PWB have not been observed in either Lake Erie or the Detroit River in more than twenty years and are believed to be extirpated (COSEWIC 2021). The Committee on the Status of Endangered Wildlife in Canada (COSEWIC) has assessed populations of PWB in Canada as Threatened due to the loss of populations in Lake Erie and the Detroit River and a continued decline in habitat quality where it remains extant (COSEWIC 2021). The Species at Risk Act (SARA) mandates the development of strategies for the protection and recovery of species that are at risk of extinction or extirpation from Canada. In response, Fisheries and Oceans Canada (DFO) has developed the recovery potential assessment (RPA; DFO 2007a, b) as a means of providing information and scientific advice. There are three components to each RPA - an assessment of species status, the scope for recovery, and scenarios for mitigation and alternatives to activities, that are further broken down into 22 elements. This report contributes to the RPA through the use of population modelling to assess the impact of anthropogenic harm to populations, identify recovery targets, and project population recovery with associated uncertainties. This work is based on a demographic approach developed by (Vélez-Espino and Koops 2009, 2012; Vélez-Espino et al. 2010). METHODS Information on vital rates was compiled to build a matrix population model (Caswell 2001) that incorporated parameter uncertainty, environmental stochasticity and density-dependence. The impact of anthropogenic harm to populations was quantified with use of elasticity and simulation analyses. Estimates of recovery targets for abundance and habitat were made with estimation of the minimum viable population (MVP) and the minimum area for population viability (MAPV). Finally, simulation analysis was used to project population abundances and make estimates of potential recovery time-frames. All analyses and simulations were conducted using the statistical program R 4.1.2 (R Core Team 2021). SOURCES DFO conducts a long-term monitoring program for Ontario unionid populations – Unionid Monitoring and Biodiversity Observation (UMBO; Sheldon et al. 2020). The monitoring program includes index sites in the Sydenham and Thames rivers. In the Sydenham River, 12 index sites were sampled between 1999 and 2015 over two sampling periods (1999–2003 and 2012– 2015). In the Thames River, 14 sites were sampled between 2004 and 2018 over two sampling periods (2004–2010 and 2015–2018). In the Sydenham River, 10 sites were sampled twice, once in each sampling period, and in the Thames River, 8 sites were sampled twice within the expected distribution of PWB. At each site, quadrat sampling was conducted following methods modified from Metcalfe-Smith et al. (2007). Sites were divided into approximately 25 blocks measuring 3 m by 5 m (15 m2). Three 1 m2 quadrats were random sampled from each block. Mussels were counted, identified to species, and measured (maximum length). In addition, shells from dead mussels were collected and aged. In total, 3,275 PWB were sampled, 3,085 in the Sydenham River and 190 in the Thames River. In addition, 61 PWB shells were aged across both rivers (35 from the Sydenham River, 25 from the Thames River, and 1 uncertain). 2 The monitoring data were analysed by van der Lee et al. (in prep.1) providing information on population trajectory, density, growth and survival rate which was incorporated into the population model. The analysis of the quadrat data was conducted using a hierarchical Bayesian model, with quadrate count modelled as a function of sample year (as a continuous variable), to estimate density, population trajectory and project population size (van der Lee et al. in prep.1). The analysis was limited to sites that were sampled twice representing years 1999–2015 in the Sydenham River and 2004–2017 in the Thames River. PWB density was greater in the Sydenham River than in the Thames River. Expected median density, projected from the fitted model, was estimated for the final year of sampling in each river. In the Sydenham River expected density in 2015 was 1.82 mussels∙m-2 (CI: 0.94–3.87) and in the Thames River expected density in 2017 was 0.12 mussels∙m-2 (CI: 0.03–0.42). Population growth, however was greater in the Thames River. PWB populations in both rivers experienced positive population growth over the respective periods (1999–2015 for the Sydenham River and 2004– 2017 for the Thames River). Population growth rates (λ) were 1.047 (CI: 1.037–1.058) in the Sydenham River and 1.157 (CI: 1.10–1.221) in the Thames River. The temporal trend was also modelled using sample period as a categorical variable, rather than year as a continuous variable, which confirmed the importance of the change in density through time. Population abundance was projected only to the survey locations (i.e., not extrapolated to other areas of the rivers). The sampled portion of the rivers covers an area of 3,600 m2 in the Sydenham River and 3,000 m2 in the Thames River. The estimated number of PWB in the sampled region of the Sydenham River was 10,504 (CI: 9,563–11,505) in 2015 and in the Thames River was 872 (CI: 696–1,091) for 2017. Additional information on Unionid life-history was taken from (Haag and Staton 2003, Haag and Rypel 2011, Haag 2012). THE POPULATION MODEL The PWB life cycle was modelled using a female only, density-dependent, birth-pulse, pre- breeding, stage-structured population matrix model with an annual projection interval (Caswell 2001, Figure 1). Figure 1. Generalized life cycle used to model the population dynamics of PWB. Fi represents stage- specific annual fertility, Pi represents the probability of surviving and remaining in stage i, and Gi represents the probability of surviving and moving to stage i+1 each year. 1 van der Lee, A.S., Goguen, M.N., McNichols-O’Rourke, K.A., Morris, T.J., and Koops, M.A. In prep. Evaluating the status and biology of an imperilled freshwater mussel, Purple Wartyback (Cyclonaias tuberculata), in Southern Ontario. In preparation. 3 The matrix consisted of 3 stages (Figure 1) representing juveniles, young adults and old adults. The juvenile stage represents immature PWB, from age-1 to the age-of-maturity (Tmat). The adult stage was divided into young and old based on expected annual somatic growth. The young adult stage represented mature individuals that continued to put resources toward somatic growth; the old adult stage represented individuals that express minimal somatic growth. The amount of growth, in mm, was estimated from the von Bertalanffy growth function (VBGF, see Equation 7 below). The division was set to the age where annual growth was < 0.5 mm (age-34). The oldest PWB identified in Canada was 92 years old (van der Lee et al. in prep.1). Therefore, stage-1 includes ages 1 to Tmat which represented 5 to 9 age classes depending on the values of Tmat (see below); stage-2 includes ages Tmat to 34 which represents 24 to 28 age-classes; and stage-3 includes ages > 34 which represents at least 58 age-classes. The projection matrix A is the product of the transition matrix B, consisting of the life-history characteristics, and the density-dependence matrix D (see Equation 10 below) representing the density-dependence effects, where: 𝐁𝐁 = � 𝑃𝑃𝑗𝑗 𝐹𝐹𝑎𝑎1 𝐹𝐹𝑎𝑎2 𝐺𝐺𝑗𝑗 𝑃𝑃𝑎𝑎1 0 0 𝐺𝐺𝑎𝑎1 𝑃𝑃𝑎𝑎2 �, (1) and: 𝐀𝐀 = 𝐁𝐁 ∘ 𝐃𝐃, (2) where the symbol ∘ represents the Hadamard product or the element by element multiplication of the matrices. Stage-specific abundance, 𝒏𝒏, each year, 𝑦𝑦, is calculated from: 𝒏𝒏𝑦𝑦+1 = 𝐀𝐀𝑛𝑛,𝑦𝑦𝒏𝒏𝑦𝑦, (3) where 𝒏𝒏 is a vector of stage-specific abundance and the population projection matrix A varies among years due to environmental conditions and population size impacting vital rates. Stage-based matrix models incorporate estimates of Fi, stage-specific fertility, Pi, the probability of survival and remaining in stage i, and Gi, the probability of surviving and moving to the next stage. Fertility, Fi, represents the number of female offspring produced per adult female annually and includes all reproductive parameters including: fecundity, sex ratio, spawning periodicity and survival from the egg state to age-1. Pi and Gi are both a function of stage-specific survival (σi) and stage-specific transition probabilities (𝜏𝜏𝑖𝑖) describing the likelihood of moving from stage 𝑖𝑖 to 𝑖𝑖 + 1, where: 𝑃𝑃𝑖𝑖 = 𝜎𝜎𝑖𝑖(1 − 𝜏𝜏𝑖𝑖) and (4) 𝐺𝐺𝑖𝑖 = 𝜎𝜎𝑖𝑖𝜏𝜏𝑖𝑖. (5) The stage transition probability (𝜏𝜏𝑖𝑖) describes the proportion of individuals in stage 𝑖𝑖 that will move to stage 𝑖𝑖 + 1, which can be estimated from (Caswell 2001): 𝜏𝜏𝑖𝑖 = 𝜎𝜎𝑖𝑖𝑇𝑇𝑖𝑖−𝜎𝜎𝑖𝑖(𝑇𝑇𝑖𝑖−1) (𝜎𝜎𝑖𝑖𝑇𝑇𝑖𝑖−1) . (6) Where 𝑇𝑇𝑖𝑖 represents the stage duration in years. Equation 6 describes the oldest individuals in stage 𝑖𝑖 moving to stage 𝑖𝑖 + 1 the following year. PARAMETERIZATION The understanding of many important characteristics of the life-history and population ecology of PWB is incomplete. When possible, model parameters were estimated directly from field data 4 (Table 1). When this was not possible parameters were estimated from generic relationships or solved for to provide a certain population state (e.g., population growth rate (λ) = 1). When a parameter was not directly estimated from PWB population specific data uncertainty was incorporated into model runs by drawing the unknown parameter from a specified probability distribution providing a range of plausible values (Table 1). This allows different model runs to incorporate different parameter combinations and provides results that span the full extent of the uncertainty. Life-History Growth of PWB was described with the von Bertalanffy growth function fitted to aged shells found in the Sydenham and Thames rivers. Ages ranged from 1 to 92 years. Length-at-age was described by (van der Lee et al. in prep.1): 𝐿𝐿𝑡𝑡 = 110.9(1 − 𝑒𝑒−0.091𝑡𝑡). (7) The growth relationship was only used to identify the age division between young and old adults which was defined as the age where growth was < 0.5 mm per year, occurring at age-34. Table 1. Parameter values incorporated in the matrix population model for PWB. Fixed parameters were held constant across model runs and uncertain parameters were sampled randomly from the specified distribution. Parameter Definition Value Fixed parameters Ma Adult instantaneous mortality rate 0.0508 σa Adult survival rate (𝑒𝑒−𝑀𝑀𝑎𝑎) 0.9504 CVM Coefficient of variation for M. Used to describe inter- annual variation in survival. 0.15 L∞ Asymptotic length (mm) 112.6 k Growth coefficient 0.0871 Ta Age division between young adult and old stages 34 ζ Generation time 25.4 years (calculated) Uncertain parameters Tmat Age-at-50%-maturity 𝑈𝑈𝑈𝑈𝑖𝑖𝑈𝑈(6, 10) σj Juvenile survival rate 𝐵𝐵𝑒𝑒𝐵𝐵𝐵𝐵(15.1, 2.9) F Fertility Solved for (Figure 3) CVF Coefficient of variation for F. Used to describe inter- annual variation in fertility. 𝑈𝑈𝑈𝑈𝑖𝑖𝑈𝑈(0.15, 0.5) α Relative fertility between young and old adults 𝑈𝑈𝑈𝑈𝑖𝑖𝑈𝑈(0.5, 2.0) λmax Maximum population growth rate 𝑈𝑈𝑈𝑈𝑖𝑖𝑈𝑈(1.1, 1.4) The age data were also used to estimate adult survival rate. Catch curve analysis was conducted using the Chapman-Robson method (Smith et al. 2012). Ages of mussels incorporated in the analysis were 7 to 92 years. Adult survival rate and instantaneous mortality estimates were 0.9504 (SE = 0.069) and 0.0508 (SE = 0.008) respectively. Haag (2012) provides a relationship to estimate instantaneous mortality for unionid species from longevity (Tmax) based on a regression analysis with data from 14 species (15 populations). The fitted relationship was 𝑀𝑀 = 4.171𝑇𝑇𝑚𝑚𝑎𝑎𝑚𝑚−1.070 (𝑟𝑟2 = 0.93). Using the maximum observed age of 92 to represent longevity the predictive relationship gives an estimate for instantaneous mortality of 5 0.035 or a survival rate of 0.966. The predicted estimate corresponds well to the catch curve analysis. The catch curve analysis provides an estimate of mean adult survival which was held constant for both young and old adult stages. No data, however, were available to inform PWB juvenile survival. To estimate juvenile survival, a relationship was fit (beta regression) between adult (𝜎𝜎𝑎𝑎) and juvenile (𝜎𝜎𝑗𝑗) survival rates based on data from 11 species reported by Haag (2012). The beta regression was fit using the betareg package in R (Cribari-Neta and Zeileis 2010). The fitted relationship was (Figure 2): 𝑙𝑙𝑙𝑙𝑙𝑙𝑖𝑖𝐵𝐵�𝜎𝜎𝑗𝑗� = 4.38𝜎𝜎𝑎𝑎 − 2.52 (𝜑𝜑 = 8.51; 𝑝𝑝 < 0.01; 𝑝𝑝𝑝𝑝𝑒𝑒𝑝𝑝𝑝𝑝𝑙𝑙 𝑟𝑟2 = 0.493). (8) Equation 8 gives a mean estimate of juvenile survival for PWB of 0.838 (variance = 0.0143). Because 𝜎𝜎𝑗𝑗 is an unknown parameter for PWB it was included in the population model as a stochastic parameter. Therefore for each model run, 𝜎𝜎𝑗𝑗 was drawn from a probability distribution to give a random value for that simulation. The distribution used to generate stochastic juvenile survival rates was 𝐵𝐵𝑒𝑒𝐵𝐵𝐵𝐵(15.1, 2.9) (Figure 2) based on the mean and half of the variance produced from Equation 8. The variance was halved to limit the potential juvenile survival rates to a reasonable range. Figure 2. Juvenile survival rate estimates used in the PWB population model. The left panel is the beta regression model fit predicting juvenile survival from adult survival with data from Haag (2012); the black line is the mean trend and the grey area are the 80% CIs. The right panel is the beta distribution used to randomly draw juvenile survival estimates in model runs (𝐵𝐵𝑒𝑒𝐵𝐵𝐵𝐵(15.1, 2.9)). The red lines indicate the mean estimate for PWB. Age-at-maturity (Tmat) was also an uncertain parameter. Jirka and Neves (1992) reported the youngest mature female PWB observed in New River, Virginia was age-6 and 58.6 mm in length. No other information on PWB maturity was available. Haag (2012) provides a relationship to estimate age-at-first-maturity from the VBGF growth coefficient (k) based on a regression analysis of 16 unionid species. The fitted relationship was 𝑇𝑇𝑚𝑚𝑎𝑎𝑡𝑡 = 0.69𝑘𝑘−1.031 − 1 (𝑟𝑟2 = 0.94). This gives an estimate for PWB, where k was 0.091, of 7.2 years. These estimates represent the age-of-first-maturity, while the model requires the value of age-of-50%- maturity. Based on the maturity ogives of other unionid species (Haag and Staton 2003) it is likely to take at least 2 years to reach 50% maturity, suggesting values in the range of 8–10 6 years. Age-at-maturity was included in the population model as a stochastic parameter that varied between model runs. Uncertainty in Tmat was represented using a discrete uniform distribution with limits of ages 6 and 10 (𝑈𝑈𝑈𝑈𝑖𝑖𝑈𝑈(6, 10)). This allowed for an equal probability for Tmat to be 6, 7, 8, 9 or 10 across model runs. PWB reproduction is largely unstudied in Canada but is similar to other members of the Unionidae family (COSEWIC 2021). Spawning likely occurs in the spring with males releasing sperm into the water and females filtering it out of the water. PWB are short-term brooders and release their glochidia (larvae mussels) in the same summer. Like all other unionids, PWB are obligate parasites of vertebrate (usually fish) hosts. Glochidia must attach to a specific host species to complete their development. The host species for PWB are Black Bullhead (Ameiurus melas), Yellow Bullhead (Ameiurus natalis), Channel Catfish (Ictalurus punctatus), and Flathead Catfish (Pylodictus olivaris) (Hove 1997). PWB attract hosts with both mantle display and an amorphous conglutinates to increase the likelihood of host infestation (Sietman et al. 2012). Mortality, however, during this phase is expected to be high with survival from the glochidial stage to settlement for most unionids in the order of 10-5 to 10-6 (Haag 2012). After the parasitic phase the juvenile mussels release from the host and settle on the river bottom. The matrix model requires an estimate of annual fertility per female. Fertility represents the number of female offspring surviving to age-1 and includes all aspects of reproduction including: fecundity, spawning periodicity, sex ratio at birth, and survival from egg to age-1. None of these aspects are known for PWB. Instead the value of fertility that would result in a stable population (population growth rate (λ) = 1) was solved for, for a given combination of life-history parameters. This provides an estimate of the average number of newly spawned females that is needed to survive to age-1 each year for the population to be stable; however, the value of each of the components of fertility remain unknown. Figure 3. Distribution of estimated fertility values (the number of female offspring surviving to age-1 produced annually per adult female) for young and old adult PWB that gave populations growth rates of 1. PWB is primarily dioecious, with little incidence of hermaphrodism and an even sex ratio (Haggerty et al. 1995). Spawning periodicity has not been measured for PWB, however, evidence from other unionid species suggests spawning may occur annually, as a high percentage (> 90%) of mature females were gravid across 6 unionid species (Haag and Staton 2003). Fecundity has not been measured for PWB. Haag (2012) provides a relationship to estimate fecundity from adult length where 𝐹𝐹𝑒𝑒𝐹𝐹𝑝𝑝𝑈𝑈𝑝𝑝𝑖𝑖𝐵𝐵𝑦𝑦 = 0.213𝐿𝐿𝑎𝑎3.146 (𝑟𝑟2 = 0.751,𝑈𝑈 = 71); 7 however, 5 outlying species were excluded from the analysis. Using an average adult size of 75 mm estimated from the geometric mean of Lmat and Lmax, calculated from the VBGF (Equation 7), the relationship predicts annual fecundity for PWB of ~170,000 eggs. Across individuals within a species, fecundity generally increases with mussel length, but there is some evidence to suggest it may plateau or even decrease at older ages, potentially indicative of reproductive senescence (Haag and Staton 2003). Therefore, it is unclear how fertility will relate between young and old adults in the population matrix model. To allow for potential differences between the two stages an additional parameter, α, was included to represent relative fertility between young and old adults. Relative fertility was incorporated as a stochastic parameter (Table 1) and drawn from a uniform distribution with limits of 0.5 and 2 ((𝑈𝑈𝑈𝑈𝑖𝑖𝑈𝑈(0.5, 2.0)). A value of 1 for α means fertility is constant throughout the adult stage. A value of α > 1 means that fertility of old adults is greater than young adults; possibly from greater fecundity due to increase body size or greater survival of glochidia from older individuals. A value of α < 1 means fertility decreases with age due to some degree of reproductive senescence. The range of α values included in the analyses allows fertility of old adults to be as little as half that of young adults (α = 0.5) and up to as high as double that of young adults (α = 2). First year survival rate is unknown. By solving for fertility, the uncertainty in the various reproductive parameters is reduced to a single value for each matrix. Across replicates with different stochastic life-history parameters the median fertility value for young adults was 0.14 (95% CI: 0.06–1.33) and for old adults was 0.17 (95% CI: 0.05–1.62). Representing the number of female offspring per female that survive to age-1 required to result in a stable population size. If fecundity is taken to be 170,000, the sex ratio is 1:1, and 95% of females spawn annually this gives an approximate survival rate from egg to age-1 of 1.9 × 10−6. A typical unionid species produces between 0.1–1.3 juveniles per year or 0.05–0.75 juvenile females per year (Haag 2012), which is in line with the estimates for PWB. Density-Dependence Density-dependence was assumed to act during the first year of life. As a largely sedentary species density-dependence is less likely to impact later life stages. For unionids, density- dependence has the potential to act at two stages during the first year of life: host infestation and settlement. PWB can use multiple host species (Channel Catfish, Black and Yellow bullhead) to complete their development which are common throughout southern Ontario (COSEWIC 2021). As a result, PWB is less limited by host availability than other unionid species (Daniel et al. 2018). In the model, it is therefore assumed that PWB are not limited by host availability. Density-dependence is assumed to act due to habitat availability during settlement following the parasitic phase. Density-dependence was implemented using a Beverton-Holt function with fertility a function of adult density: 𝑝𝑝 = 𝐹𝐹𝑚𝑚𝑎𝑎𝑚𝑚 𝐹𝐹1� 1+𝑏𝑏 𝐾𝐾� 𝑁𝑁𝑎𝑎 , (9) where Na is adult density, K is carrying capacity, b is the density-dependence coefficient, F1 is fertility when the population is stable and Fmax is fertility when the population is at maximum population growth rate (λmax). Maximum population growth rate was unknown and included as a stochastic parameter (Table 1). λmax was drawn from a uniform distribution with a minimum value of 1.1 and maximum value of 1.4 (𝑈𝑈𝑈𝑈𝑖𝑖𝑈𝑈(1.1, 1.4)). The maximum value was based on the greatest observed P/B ratio reported for a unionid species (Patterson 1985). Fmax was then solved for in the same manner as fertility to find the value that results in a population growing at its maximum rate. The value b was calculated by rearranging Equation 9, as: 𝑏𝑏 = 𝐹𝐹𝑚𝑚𝑎𝑎𝑚𝑚 𝐹𝐹1 − 1⁄ . The density-dependence matrix, D, was structured as: 8 𝐃𝐃 = � 1 𝑝𝑝 𝑝𝑝 1 1 1 1 1 1 �, (10) and incorporated into Equation 2. Stochasticity Stochasticity was incorporated into simulations at two levels: to account for parameter uncertainty between model runs and to account for environmental variability within model runs. Parameter uncertainty was implemented by drawing uncertain parameters from defined probability distributions prior to each model run (Table 1). Uncertain parameters included: juvenile survival rate, age-at-maturity, relative fertility, maximum population growth rate, and the inter-annual variance in fertility. Environmental stochasticity was incorporated into the model by allowing certain parameters to vary between years within a model run to simulate the changes in vital rates that occur naturally due to variation in environmental conditions. Environmental stochasticity was applied to annual survival rate and fertility. Juvenile and adult survival were allowed to vary independently following a log-normal distribution with a CV of 0.15, which was the approximate amount of error in annual survival rate estimated from the Chapman-Robson catch curve analysis. The variability in fertility was unknown. Therefore, the amount of variability in fertility incorporated into simulations was included as a stochastic variable (Table 1). Fertility was varied using a log- normal distribution with a CV drawn from a uniform distribution with limits of 0.15 and 0.5 (𝑈𝑈𝑈𝑈𝑖𝑖𝑈𝑈(0.15, 0.5)). The minimum value, 0.15 was chosen to mirror the variability in survival rate as first year survival is an element of fertility. The upper limit was chosen arbitrarily but is set to allow greater variation in recruitment success between years. At the upper 95% confidence limit of old adult fertility, 1.62, a CV of 0.5 allows for an upper 95% confidence limit value for fertility of 3.66 meaning that under exceptional environmental conditions 3.66 female offspring could survive to age-1 per adult female in populations with high juvenile mortality and later maturity (conditions that require greater fertility for population stability). IMPACT OF HARM The impact of anthropogenic harm to a PWB population was assessed with deterministic elasticity analysis of the projection matrices and stochastic simulations. Elasticity analysis of matrix elements provides a method to quantify the impact of changes to vital rates on a population. Specifically, elasticities measure the proportional change to population growth rate (λ) that results from a proportional change in a vital rate (v). For example, an elasticity of λ value of 0.2 for juvenile survival indicates that a 10% change in juvenile survival rate (e.g., 0.84 × (1 + 0.1) = 0.92) would cause a 2% increase in population growth rate (e.g., 1 × (1 + 0.1 × 0.2) = 1.02). Elasticities are useful as they allow for assessment of how impactful changes to vital rates and other model parameters are to a population and because they represent proportional changes their values are directly comparable. They are preferable to simulation analyses because of the speed with which they can be estimated allowing for many more perturbations to be examined than simulations. Elasticities are limited, however, as they represent permanent changes, accurately represent only small perturbations (i.e., < 30% changes), and assume all other model parameters remain unchanged. Therefore, simulation analysis was used to examine the effects of transient or periodic harm to a population. 9 Elasticity of λ Elasticities of λ (ελ) are calculated by taking the scaled partial derivatives of λ with respect to a vital rate (𝜈𝜈, Caswell 2001): 𝜀𝜀𝜆𝜆 = 𝜈𝜈 𝜆𝜆 ∑ 𝜕𝜕𝜆𝜆 𝜕𝜕𝑎𝑎𝑖𝑖,𝑗𝑗 𝜕𝜕𝑎𝑎𝑖𝑖,𝑗𝑗 𝜕𝜕𝜈𝜈𝑖𝑖,𝑗𝑗 , (11) where aij is the projection matrix element in row i and column j. A range of potential vital rate elasticities were estimated by allowing uncertain parameters (Table 1) to vary. As well, population growth rate was allowed to vary by incorporating the effect of density on fertility (Equations 9 & 10). This was accomplished by randomly generating adult density as a proportion of carrying capacity (0–1000%) using a draw from a uniform distribution on a log scale (𝑒𝑒𝑈𝑈𝑛𝑛𝑖𝑖𝑈𝑈(log (0.01),log(𝐾𝐾×10))). λ = 1 when Na = K, λ < 1 when Na > K and λ > 1 when Na < K up to a maximum λ of 1.4. The elasticity analysis was replicated 5,000 times. The maximum allowable harm can be estimated from elasticities by calculating the change in a vital rate that allows a population to maintain a population growth rate ≥ 1. Allowable harm applies when a population has an initial λ > 1. Maximum allowable harm is estimated as (Vélez- Espino and Koops 2009): 𝑀𝑀𝐵𝐵𝑀𝑀𝑖𝑖𝑀𝑀𝑝𝑝𝑀𝑀 𝐵𝐵𝑙𝑙𝑙𝑙𝑙𝑙𝑎𝑎𝐵𝐵𝑏𝑏𝑙𝑙𝑒𝑒 ℎ𝐵𝐵𝑟𝑟𝑀𝑀 = � 1 𝜀𝜀𝜆𝜆 � �1−𝜆𝜆 𝜆𝜆 �. (12) Maximum allowable harm is estimated for the Sydenham and Thames rivers based on the full range of estimates of their population growth rates (van der Lee et al. in prep.1); Sydenham River: 1.03–1.07, Thames River: 1.07–1.27. Simulation Simulation analysis was used to investigate the impacts of periodic stage-specific harm on adult population density. Life-stages were impacted by a level of harm, ranging from 0 to 99%, at different frequencies: 1, 2, 5, and 10 years, over a 100 year simulation. The initial carrying capacity was then compared to the mean population size over the final 15 years of the simulation to determine the effect of the harm, quantified as the proportion of initial Ka. The frequency indicates how often harm was applied. A frequency of 1 indicates that harm is constant and applied every year, where a frequency of 10 indicates that harm is periodic and applied once every 10 years. As a density-dependent model it is assumed that the population is able to recover in between applications of harm as conditions are returned to the initial state and no competitors exist as this is a single-species model. The simulations incorporated parameter uncertainty and environmental stochasticity and were replicated 1,000 times. RECOVERY TARGETS Abundance: Minimum Viable Population (MVP) The concept of demographic sustainability was used to identify potential minimum recovery targets for PWB. Demographic sustainability is related to the concept of a minimum viable population (MVP, Shaffer 1981), and was defined as the minimum adult population size that results in a desired probability of persistence over 250 years (~ 10 PWB generations, where generation time (ζ) was estimated from stochastic projection matrices (Caswell 2001) with λ = 1 where ζ = 25.4 (CI: 21.9–29.1)). MVP was estimated using simulation analysis which incorporated, parameter uncertainty, environmental stochasticity and density-dependence. 10 Important elements incorporated in population viability analysis include: the choice of time frame over which persistence is determined, the severity and frequency of catastrophic events, and the quasi-extinction threshold below which a population is deemed unviable. The choice of time frame is arbitrary and without biological rational. It must be long enough to be of consequence for the population examined but not so long as to be an unreasonable timeframe for management considerations. 250 years was selected as it equates to approximately 10 PWB generations. Catastrophes represent an event that causes a > 50% decrease in population size. A PWB population may be impacted by catastrophes such as severe drought or disease (Haag 2019). The rate that severe catastrophes impact PWB populations is unknown. Based on a meta- analysis, Reed et al. (2003) determined that among vertebrate populations, catastrophic die-offs (defined as a one-year decrease in population size > 50%) occurred at a rate of 14%∙generation-1 on average. The rate of catastrophic die-offs among invertebrate populations is unknown and was included in simulations as a stochastic parameter by drawing from a uniform distribution with limits of 0.05 and 0.2 (𝑈𝑈𝑈𝑈𝑖𝑖𝑈𝑈(0.05, 0.2)) to allow for catastrophe rates of 5% to 20% per generation. This equates to a large scale die-off occurring approximately every 507 years (at a 5%∙generation-1 rate) to every 126 years (at a 20%∙generation-1 rate), calculated by dividing generation time by the catastrophe rate. The impact of a catastrophe affects all life- stages simultaneously and was drawn randomly from a beta distribution scaled between 0.5 and 1 with shape parameters of 0.762 and 1.5 (Reed et al. 2003, 𝐵𝐵𝑒𝑒𝐵𝐵𝐵𝐵(0.762, 1.5) × (1 − 0.5) + 0.5), representing the probability of a 50 to 100% decline in population size. Quasi-extinction represents the compounding effects of Allee effects, demographic stochasticity and inbreeding depression (Lande 1988) leading a population to extinction once the threshold is crossed. The value of the quasi-extinction threshold cannot be empirically measured; therefore, 25 adult females was used as a reasonable approximation (Morris and Doak 2002). Simulations were conducted for populations with various initial abundances of adult females, ranging from 50 to 5,000 (initial density represented carrying capacity, K, where λ = 1) and catastrophe rates ranging from 5% to 20% per generation. The simulations were conducted for 250 years and replicated 10,000 times. The number of quasi-extinctions was tracked across simulations and the probability of extinction (Pext) was modelled with a logistic regression that was a function of log10 transformed initial female density (log10(Na)), catastrophe rate (Pcat), and all uncertain population model parameters (Table 1) where: 𝑃𝑃𝑒𝑒𝑚𝑚𝑡𝑡 = 1 1+𝑒𝑒−(𝑿𝑿𝛽𝛽), (13) Where X is a matrix of all covariates and β a vector of coefficients including the intercept. Habitat: Minimum Area for Population Viability (MAPV) Minimum area for population viability (MAPV) represents the quantity of habitat required to support a population of MVP size. MAPV is estimated by dividing the MVP estimate by density. Density estimates were available for PWB populations in the Sydenham and Thames rivers from a hierarchical Bayesian model fit to quadrat survey data. Median density, estimated from the model fit, in the Sydenham River was 1.82 mussels m-2 (CI: 0.94–3.87) in 2015, and in the Thames River was 0.12 mussels m-2 (CI: 0.03–0.42) in 2017 (van der Lee et al. in prep.1). These data represent whole population density, whereas MVP is specific to adult females. A 1:1 sex ratio can be assumed and based on the length frequency distribution (van der Lee et al in prep.1) sampled mussels were approximately 87% adult in the Sydenham River and 55% adult in the Thames River. In addition, both these populations have a positive population trajectory with population growth rate estimates of 1.047 (CI: 1.037–1.058) in the Sydenham River and 11 1.157 (CI: 1.10–1.221) in the Thames River. Therefore, the current densities do not represent carrying capacity for the habitat and would therefore produce over-estimates of the quantity of habitat required to house an MVP sized population, which assumes population stability. The matrix population model can be used to provide an estimate of the density what would give a stable population size based on the current density and population growth rate. The density that gives a stable population size was determined and used to provide estimates of MAPV. RECOVERY TIMES Recovery time simulations were conducted for Thames River populations. Current population size was taken from the quadrat model projection and represents the number of PWB in the sampled portion of the river in 2017 (872, CI: 696–1,091, 3,000 m2, van der Lee et al. in prep.1). This value was converted to the number of adult females assuming a 1:1 sex ratio and 55% adult sample. Current population trajectory was also taken from the quadrat model (λ = 1.157, CI: 1.10–1.221). Carrying capacity of the available habitat was solved for using the projection matrix given the current abundance and λ value based on the assumed density-dependence relationship. The time taken for the population to exceed MVP for a given catastrophe rate was recorded. Simulations were replicated 10,000 times. RESULTS IMPACT OF HARM The impact of anthropogenic harm to PWB populations was assessed with two analyses: deterministic elasticities and simulations. Elasticity of λ Elasticities demonstrate which vital rates, if changed, have the greatest impact on a population’s growth rate. In most instances, PWB populations were most sensitive to changes to adult survival rate (Figure 4). Juvenile survival rate had the next greatest elasticities and under some conditions had a greater elasticity than adult survival. Age-at-maturity produced negative elasticities indicating that increases to Tmat would cause a decrease in λ. Fertility encompasses all aspects of reproduction, including egg production and first year survival rate, and the elasticity value for fertility can apply to any of these aspects independently. Under most conditions, however, fertility had small elasticities relative to adult survival. As well, the elasticity of fertility of old adults was small and less than that of young adults. 12 Figure 4. Elasticity of λ analysis results for PWB populations represented with violin and box plots. Results reflect uncertainty in life-history characteristics and different values of population growth. F represents fertility, σ represents survival rate (juvenile, j, and adult, a), and Tmat represents age-at- maturity. Table 2. Pearson correlations between uncertain life-history parameters and elasticity estimates. Absolute values > 0.1 indicate when elasticity estimates were sensitive to a the value of a life-history parameter. Life-History Parameter Vital Rate Elasticity Lambda Age-at-maturity Relative Fertility Juvenile Survival Juvenile Survival 0.91 0.26 -0.06 -0.14 Adult Survival -0.91 -0.07 0.30 -0.12 Fertility – Young Adult 0.99 -0.05 -0.10 -0.08 Fertility – old adult -0.03 0.19 0.76 -0.10 Fertility - Total 0.99 -0.03 -0.06 -0.08 Age-at-maturity -0.78 -0.11 0.05 0.52 The elasticity values for vital rates were influenced by population state (λ) and other life-history parameter values (Figure 5). Correlations between elasticity estimates and life-history parameter values indicated how influential the life-history parameter values are on the elasticity estimates (Table 2). Absolute correlation values > 0.1 indicate that the elasticity estimate was sensitive to the value of a particular life-history parameter. The most influential parameters were population growth rate, age-at-maturity, and relative fertility (Table 2). The elasticity values of young adult fertility and juvenile survival increased significantly with λ while that of adult survival and age-at-maturity decreased; absolute correlation values were > 0.9 indicated a very strong dependency. When λ > ~1.2 the elasticity of juvenile survival was greater than adult survival, while the adult survival elasticity was greater when λ < ~1.075. Between those values other life history parameters (especially age-at-maturity and relative fertility) had influence over which elasticity was greater. The elasticity of juvenile survival was also influenced by age-at-maturity. When maturity occurred later, changes to juvenile survival rate had a greater influence on λ and at greater λ values the importance of Tmat increased. The elasticity value of adult survival was influenced by relative fertility. With increased relative fertility, a greater proportion of reproduction occurs later in life, which caused an increase in population sensitivity to adult 13 survival. Relative fertility influenced the elasticity value of old adult fertility as well, however, the magnitude of the difference was small (~0.01). The value of juvenile survival only had an influence on the elasticity values of age-at-maturity. The maximum amount of harm consistent with the objective to maintain stable or growing populations (i.e., population growth rate at or above 1) was estimated for the Sydenham and Thames river systems based on their current states (Table 3). The PWB populations in the two rivers differ in their current rates of population growth and, therefore, the estimates of maximum allowable harm differ. As well, this leads to different patterns in elasticity estimates such that in the Sydenham River the adult survival produced the lowest stage-specific maximum allowable harm estimate while in the Thames River juvenile survival produced the lowest stage-specific maximum allowable harm estimate. Neither system was particularly sensitive to fertility with lower confidence interval estimates of maximum allowable harm > 49 and 81% for the Sydenham and Thames Rivers respectively. Maximum allowable harm if both juvenile and adult stages were impacted was also estimated (Table 3). Figure 5. Elasticity of λ analysis results for PWB populations plotted against population growth rate. Panels represent the vital rate elasticities. Colour indicates the relative fertility parameter, α, and plot symbol represents age-at-maturity. The vertical lines represent mean estimate of population growth rate for Sydenham River (1.05, solid) and Thames River (1.16, dashed) populations. NOTE: different y-axis scales. 14 Simulation Simulation analysis was used to investigate the impact of anthropogenic harm to stable population size and investigate the effects of periodic perturbations occurring annually (for comparison to elasticity analysis), every second year, fifth year, and tenth year (Figure 6). Harm applied to the adult stage or all stages had a significant impact on abundance. Small amounts of mortality had large impacts on population size. Annual mortality rates of only 1–2% applied to the adult stage resulted in a 25% reduction in population size. With biennial harm this became ~3%, ~7% with harm occurring every 5 years, and ~15% every 10 years. The effects of harm were smaller but still significant if only juvenile mussels were impacted. PWB populations were not as impacted by harm to fertility, representing interruptions to reproduction or harm to glochidia/post-settlement age-0 mussels, particularly if harm occurred infrequently. Table 3. Maximum allowable harm estimates for Sydenham River and Thames River populations of PBW. The values represent the maximum percent decrease in vital rates that will allow the population to maintain a population growth rate ≥ 1. Allowable harm estimates are made based on estimated population growth rates for each river (van der Lee et al. in prep.1); Sydenham River: 1.03–1.07, Thames River: 1.07–1.27. LCI and UCI are the lower and upper confidence intervals respectively. Max. Allowable Harm Vital Rate Median LCI UCI Sydenham River Juvenile Survival 9.9 6.1 16.1 Adult Survival 6.2 3.7 10.0 Juvenile & Adult Survival 3.8 2.5 5.5 Fertility 65.4 49.3 84.6 Thames River Juvenile Survival 19.4 11.0 31.3 Adult Survival 24.5 9.2 48.2 Juvenile & Adult Survival 10.9 5.4 17.6 Fertility > 100 80.7 > 100 RECOVERY TARGETS Abundance: Minimum Viable Population (MVP) Demographic sustainability was assessed using simulation analysis which incorporated, parameter uncertainty, environmental stochasticity and density-dependence. Simulation outputs, binomial quasi-extinctions (1: extinct; 0: extant), were fitted using a logistic regression as a function of adult female population size, catastrophe rate and uncertain life-history parameters. Only population size and catastrophe rate had a significant effect on extinction probability (Table 4, Figure 7). None of the life-history parameters had a significant influence on extinction risk over the range of values included (Table 1) and over the length of time (250 years; ~ 10 generations) simulated. 15 Figure 6. Results of simulation analysis examining the proportional change in stable population size that resulted for anthropogenic harm occurring at different frequencies (1, 2, 5, and 10 year intervals). Harm was applied to different life stages: Fertility, juvenile, adult, and all life stages, represented by colour. The horizontal line indicates a 25% reduction in population size. Minimum viable population (MVP) size was estimated from the logistic regression model using randomly selected catastrophe rates (𝑃𝑃𝑐𝑐𝑎𝑎𝑡𝑡 = 𝑈𝑈𝑈𝑈𝑖𝑖𝑈𝑈(0.05, 0.20)) and finding the adult female population size that resulted in a 1% probability of extinction. The median MVP estimate was ~1,400 (CI: 950–2,000). If a 1:1 sex ratio is assumed then MVP including all adult PWB was ~2,800 (CI: 1,900–4,000). Extinction probability of any adult female population size can be estimated from the fitted logistic relationship for a given catastrophe rate: 𝑃𝑃𝑒𝑒𝑚𝑚𝑡𝑡 = 1 1+𝑒𝑒−(5.41−3.50 𝑙𝑙𝑙𝑙𝑙𝑙10(Na)+8.03). (14) 16 Table 4. Logistic regression model results for PWB extinction probability. Na represents adult female population size and Pcat represents catastrophe rate. Parameter Value Standard error p-value Intercept 5.41 0.44 <0.001 log10(Na) -3.50 0.15 <0.001 Pcat 8.03 1.78 <0.001 Figure 7. The probably of quasi-extinction as a function of adult female populations size. The rug plot indicates quasi-extinctions. The solid line represent the logistic regression trend with the grey region representing the confidence intervals (for Pcat ranging from 5% to 20% per generation). The dashed line indicates a 1% probability of extinction. Habitat: Minimum Area for Population Viability (MAPV) The habitat quantity required to support an MVP sized population of PWB was estimated from the expected density of adult females in a stable population. This was based on the estimate of current population density, population growth rate, and the projection matrix with the assumed density-dependence relationship. Expected density of adult females in the Sydenham River that gave a λ = 1 was 2.21 mussels∙m-2 (CI: 1.14–5.04) and for the Thames River was 0.48 mussels∙m-2 (CI: 0.08–4.47). This corresponds to MAPV estimates of 623.3 m2 (CI: 251.9– 1,396.9) and 2,900 m2 (CI: 301.5–17,166.3) for the Sydenham and Thames Rivers respectively (Figure 8). RECOVERY TIMES Time to recovery was estimated for Thames River population as the time taken for the current population to reach MVP size. Current population size applies only to a 3,000 m2 sampled portion of the river. Median time to recovery was 20 years (CI: 10-240) (Figure 9). The Sydenham River population currently exceeds the estimated MVP. 17 Figure 8. Density plot of the estimated minimum area for population viability (MAPV) for PWB populations in the Sydenham River (red) and Thames River (blue). The vertical lines indicate the median estimates. Figure 9. The distribution of recovery times for Thames River PWB populations. DISCUSSION A population model for PWB was created to make predictions on how a population may respond to anthropogenic harm, estimate potential recovery targets for abundance and habitat, and estimate recovery times for depressed populations. The model accounted for parameter uncertainty, environmental stochasticity and density-dependent effects. When the population was stable (λ = 1) the elasticity value for adult survival ranged from ~0.65– 0.92. This indicates that, for example, a 5% mortality applied to adult PWB could result in population growth rate decreasing to 0.968–0.954. Simulation analysis was used to determine 18 the effect of mortality on population size while accounting for density-dependence. With this same 5% mortality applied annually, the population decreased by an average ~45%. This is a large effect from a small level of mortality indicating a high degree of population sensitivity to potential harm to adults. When the population was growing (λ > 1) it became increasingly sensitive to juvenile mortality. When λ ≥ ~1.2 the elasticity value for juvenile survival was generally greater than the elasticity value for adult survival. Population growth for PWB in the Thames River was 1.157 (CI: 1.10–1.221) (van der Lee et al. in prep.1), therefore, harm applied to juveniles may have an equal or greater effect on population state as harm applied to adult PWB. Typically, changes in fertility did not have as large an impact on population growth rate or stable adult density as changes to juvenile or adult survival. The pattern of increased sensitivity to the adult stage and low sensitivity to fertility has been identified in other long-lived unionids (Haag 2012). Elasticity estimates were impacted by the value of uncertain life-history parameters. For example, the value of age-at-maturity impacted the elasticity estimate for juvenile survival. As Tmat increased the population generally became more sensitive to the juvenile stage. For example with Tmat = 6, mean 𝜀𝜀𝜆𝜆 for juvenile survival across population growth rates was ~0.37 while with Tmat = 10 it was ~0.63. The elasticity for adult survival generally increased with relative fertility, α; the adult stage became more important as a greater proportion of reproduction occurred later in life. As more research is conducted on PWB and the questions around certain life-history parameters are answered, more exact elasticity values can be extracted from these results. Estimates of maximum allowable harm were made for PWB populations in the Sydenham and Thames rivers based on their current estimates of population growth rate and accounting for uncertainty in life-history parameters. The estimates represent the level of mortality that would cause population growth rate to decrease to 1. Because there is uncertainty in life-history parameters and the estimates of population growth it is prudent to take the lower confidence interval as the representation of maximum allowable harm. Application of maximum allowable harm will remove all the surplus production from the population that is generating positive population growth. The consequence of this will be to stop population growth, and if the population is below a recovery target, the population will never reach that target (i.e., recovery will be jeopardized). Estimates of potential recovery targets for population abundance were provided using simulation analysis to determine the population size required for demographic stability through estimates of minimum viable population size. The median estimate of the number of adult females required for a 1% chance of extinction over 250 years was ~1,400 (CI: 950–2,000). This equates to ~2,800 (CI: 1,900–4,000) adults if a 1:1 sex ratio is assumed. The estimate of MVP was significantly influenced by the catastrophe rate included in simulations. Lower MVP estimates are representative of lower frequencies of catastrophes (5%∙generation-1) and higher MVP estimates are representative of greater frequencies of catastrophes (20%∙generation-1). Reed et al. (2003) conducted a meta-analysis to determine the frequency of catastrophes across taxa and found an average rate of 14%∙generation-1. This analysis did not include unionid species or species with generation times as long a PWB. The frequency of catastrophes is an unknown variable and was very impactful on model results. Ultimately, it represents how stable the environment is expected to be over long timeframes. The rate of catastrophe may vary among locations and therefore the most appropriate recovery target may also vary. Estimates based on more frequent catastrophes are more conservative especially if the frequency of large scale stochastic disturbances increases with climate change. Uncertain life-history characteristics did not have an impact on MVP estimates. This is likely due to the manner in which fertility was determined after selecting other life-history parameters: 19 juvenile survival, age-at-maturity, and relative fertility. Fertility was determined by finding the value that gave a stable population and therefore the same number of mature adults would result regardless of life-history characteristics. As well, the MVP simulations only represented ~10 generations (250 years) for PWB. Preliminary analysis demonstrated that with longer simulations (e.g., 500 years) some life-history characteristics, such as maximum population growth rate and age-at-maturity, can influence persistence probability and MVP size as they impact how quickly a population can respond to a perturbation. Estimates of population size were made for Sydenham River and Thames River populations (van der Lee et al. in prep.1) for the area sampled in each river. These represent 3,600 m2 and 3,000 m2 sections of each river respectively. The estimated number of PWB in the sampled region of the Sydenham River was 10,504 (CI: 9,563–11,505) for the year 2015 and in the Thames River was 872 (CI: 696–1,091) for the year 2017. Therefore the current Sydenham River population is larger than MVP after accounting for juveniles in the sample. The Thames River population estimate was < MVP; however, PWB undoubtedly inhabit other areas of the river and therefore the population, if taken as the whole river, is likely larger than the estimate. Regardless, given the estimated level of population growth the Thames River population may reach MVP size in the survey locations in 20 years (CI: 10–240). Based on the current abundance projections (converted to adult females) for the Sydenham and Thames Rivers, the extinction risk for each over 250 years was 0.16% (CI: 0.09–0.28) and 15.3% (CI: 7.8–25.7) respectively. These likely represent over-estimates of extinction risk as they assume the populations are stable and that there are no other PWB at locations outside of the surveyed area of each river. Estimates of abundances for MVP were converted into habitat requirements by dividing MVP by expected stable population densities. Population densities of PWB in the Sydenham and Thames rivers were estimated from quadrat data. Median density in the Sydenham River was 1.82 mussels m-2 (CI: 0.94–3.87) in 2015, and in the Thames River was 0.12 mussels m-2 (CI: 0.03–0.42) in 2017. These populations are growing and therefore the density estimates do not reflect carrying capacity of the habitat. The matrix model, estimates of population growth, and estimates of population density were used to estimate population density of a stable population in each system: 2.21 mussels∙m-2 (CI: 1.14–5.04) in the Sydenham River and 0.48 mussels∙m-2 (CI: 0.08–4.47) in the Thames River. These estimates are specific to adult females. The estimate for the Thames River is more uncertain because current densities are smaller and further from carrying capacity. The corresponding MAPV estimates were 623.3 m2 (CI: 251.9– 1,396.9) and 2,900 m2 (CI: 301.5–17,166.3) for the Sydenham and Thames Rivers respectively. Upper confidence interval estimates correspond to lower estimates of current density and population growth rate and greater catastrophe rates. In addition to the Sydenham and Thames Rivers, PWB inhabits the Ausable River. Many aspect of the analysis may apply to the Ausable River population if life-history characteristic, such as growth and mortality, are similar to the other southwestern Ontario populations. The estimates of elasticity and MVP could apply to the Ausable River population. An estimate of population growth rate estimated for the Ausable River population indicate the population is likely stable. If this is the case then there is little to no scope for harm as any disturbance will likely cause population decline and increase the chances of extirpation. MAPV estimates were specific to density and population growth rate estimates for the Sydenham and Thames Rivers but could be made for the Ausable river with river-specific densities. UNCERTAINTIES There were significant uncertainties in the parameterization of the population model for PWB. Uncertain parameters included: nearly all aspects of fertility, juvenile survival, age-at-maturity, 20 and maximum population growth rate. Rather than selecting specific values for these parameters that may be representative of PWB, a range of potential values were used, represented by probability distributions. This allows all the potential life-history dynamics of PWB populations to be represented in the simulations and gave a range of results. Inclusion of parametric uncertainty is important and can have significant influence on the conclusions drawn from population viability analysis (McGowan et al. 2011). As additional information about PWB life-history becomes available more precise results can be extracted from the simulations. The density-dependence relationship was assumed and may have influenced model results. Density-dependence was assumed to follow a Beverton-Holt relationship and act in the first year of life with fertility a function of adult female density. This was meant to represent settlement success following the parasitic phase of glochidia. It was assumed that, as sedentary animals, density-dependence may not act later in life, and that because PWB use multiple host species, which are relatively common, they may not be host limited. In addition, PWB sensitivity to perturbation in fertility were relatively small. This value is likely representative of the effects that perturbations to host density would have on the population. For estimates of MVP, what is important is whether or not density-dependence acts and the strength of density-dependence, more so than how it acts. For example, the value of maximum population growth rate has been found to influence MVP estimates (Lamothe et al. 2021) while the shape of the density- dependence relationship did not (van der Lee and Koops 2021). Inclusion of density- dependence allows populations to have some degree of recovery following large perturbations leading to much smaller estimates of MVP than when density-dependence is omitted (e.g., Roberts et al. 2016). PWB in each of the Sydenham and Thames rivers were treated as a single population, with panmictic reproduction and subject to the same environmental conditions. More complex population structure is possible which may impact how abundance is estimated and the persistence probability of the population as a whole. Meta-population structure, with migration or inter-breeding among sub-populations, can drastically increase population persistence but this decreases when stochastic environmental conditions among population are correlated (Palmqvist and Lundberg 1998, Reed 2004). If there is a degree of independence among PWB sub-populations within the Sydenham or Thames rivers persistence probability may be greater than estimated here. The interaction between PWB and its host species was not explicitly incorporated into the population model though likely has a significant influence on PWB population ecology. PWB may be able to parasitize multiple catfish host species which are currently common in their shared habitat (COSEWIC 2021). As a result, PWB is less limited by host availability than other unionid species (Daniel et al. 2018). However, the possibility exists that impacts to host species populations could significantly influence PWB persistence or recovery. The abundance of unionid species in a river system tends to be correlated with the abundance of their host species (Haag 2012). As well, previous simulation analysis of Hickorynut (Obovaria olivaria) populations demonstrated that host abundance and trajectory was a strong determinate of mussel abundance (Young and Koops 2013). The population dynamics of PWB may therefore be related to the population dynamics of their hosts. For example, a host population may be subject to random catastrophes in the same manner as was assumed for PWB in MVP simulations. Even if a PWB population avoids the impacts of the catastrophic mortality initially they could be greatly impacted if their hosts become limited. As well, PWB recovery would be limited by host abundance. These possibilities were not included in the simulation but could result in larger MVP estimates and longer recovery times. 21 ELEMENTS Element 3: Estimate the current or recent life-history parameters for PWB. Life-history parameters for PWB were estimated for populations in the Sydenham and Thames rivers in van der Lee et al. (in prep.1). Values applied to the population model are summarized in Table 1. Many life-history characteristics of PBW remain unknown. Uncertainties were represented with probability distributions to capture the range of values possible for the species (Table 1). The Methods section provides explanation of the parameter selections. Element 12: Propose candidate abundance and distribution target(s) for recovery Candidate abundance targets were estimated using population viability analysis and estimates of minimum viable population (MVP). Simulations incorporated density-dependence, environmental stochasticity, parameter uncertainty, and random catastrophes. Persistence probability was influenced by population size and catastrophe rate. The population size required to provide a 1% extinction probability over 250 years was ~2,800 (CI: 1,900–4,000). Larger MVP estimates correspond to simulations with more frequent catastrophes. Element 14: Provide advice on the degree to which supply of suitable habitat meets the demands of the species both at present and when the species reaches the potential recovery target(s) identified in element 12. The quantity of habitat required to support an MVP sized population of PWB was estimated by determining the density of a stable population given the population’s current density and growth rate. Estimated adult female densities of a stable population were: 2.21 mussels∙m-2 (CI: 1.14– 5.04) in the Sydenham River and 0.48 mussels∙m-2 (CI: 0.08–4.47) in the Thames River. This corresponds to MAPV estimates of 623.3 m2 (CI: 251.9–1,396.9) and 2,900 m2 (CI: 301.5– 17,166.3) for the Sydenham and Thames rivers respectively. Upper confidence interval estimates correspond to lower estimates of current density and population growth rate and greater catastrophe rates. The quantity of habitat in the Sydenham and Thames rivers available to PWB exceed the MAPV estimates. The MAPV estimate only identifies the amount of habitat required to house an MVP sized population and does not account for other considerations such as rearing habitat or spatial configuration of populations. For example, unionid mussels reproduce by spermcasting where males release sperm which is siphoned from the water by females during filter feeding (Haag 2012). In a lotic environment, successful fertilization can only be achieved when males exist upstream of females and only habitat where successful fertilization is possible can contribute to MAPV. Element 15: Assess the probability that the potential recovery target(s) can be achieved under the current rates of population dynamics, and how that probability would vary with different mortality (especially lower) and productivity (especially higher) parameters. Population count was estimated from the quadrat data model for the Sydenham and Thames rivers (van der Lee et al. in prep.1) but was restricted to the survey area. The estimated number of PWB in the sampled region of the Sydenham River was 10,504 (CI: 9,563–11,505) in 2015 and in the Thames River was 872 (CI: 696–1,091) for 2017. The estimate for the Sydenham River was greater than MVP. The estimate for the Thames River was less than MVP; however, PWB occupy other areas of the river. Simulations were conducted to determine how long it would take for the Thames River population (survey area 22 only) to reach MVP size given its current trajectory. The results indicated that the Thames River population may reach MVP size in the survey location in 20 years (CI: 10–240). Element 19: Estimate the reduction in mortality rate expected by each of the mitigation measures or alternatives in element 16 and the increase in productivity or survivorship associated with each measure in element 17. No clear links have been identified between the mitigation measures and PWB mortality rates or productivity. Element 20: Project expected population trajectory (and uncertainties) over a scientifically reasonable time frame and to the time of reaching recovery targets, given mortality rates and productivities associated with the specific measures identified for exploration in element 19. Include those that provide as high a probability of survivorship and recovery as possible for biologically realistic parameter values. Population growth rates of PWB populations in the Sydenham and Thames Rivers were estimated from quadrat survey data collected between 1999–2015 and 2004–2017 respectively (van der Lee et al. in prep.1). Population growth rates (λ) were estimated to be 1.047 (CI: 1.037– 1.058) in the Sydenham River and 1.157 (CI: 1.10–1.221) in the Thames River. Element 21: Recommend parameter values for population productivity and starting mortality rates and, where necessary, specialized features of population models that would be required to allow exploration of additional scenarios as part of the assessment of economic, social, and cultural impacts in support of the listing process. The parameter values incorporated in the population models are based on the best available data for PWB in Canada and should be used for any future population modelling. Details regarding how the parameters were estimated and source data used are outlined in the Methods section of this report. Element 22: Evaluate maximum human-induced mortality and habitat destruction that the species can sustain without jeopardizing its survival or recovery. The impact of anthropogenic harm was assessed with deterministic elasticity analysis and model simulations. PWB populations were highly sensitive to perturbations to the adult stage under most conditions. At greater population growth rates (> ~1.2), however, PWB became more sensitive to perturbations to the juvenile stage. Maximum human-induced mortality was estimated for the Sydenham River and Thames River populations based on their current rates of populations growth. The estimates represent the percent mortality that would cause population growth rate to decrease to 1 (Table 3). As there is a high degree of uncertainty in the model it is precautionary to use the lower confidence interval of the harm assessment. The population in the Ausable River low density and not growing (COSEWIC 2021), therefore there is likely little to no scope for additional anthropogenic harm. REFERENCES CITED Caswell, H. 2001. Matrix population models: construction, analysis, and interpretation. Sinauer Associates, Sunderland, MA. 722 p. COSEWIC. 2021. COSEWIC assessment and status report on the Purple Wartyback (Cyclonaias tuberculata) in Canada. Committee on the Status of Endangered Wildlife in Canada, Ottawa, ON. xi + 64 pp. Cribari-Neta, F., and Zeileis, A. 2010. Beta regression in R. J. Stat. Softw. 34(2): 1–24. https://publications.gc.ca/site/eng/9.901774/publication.html https://publications.gc.ca/site/eng/9.901774/publication.html 23 Daniel, W.M., Cooper, A.R., Badra, P.J., and Infante, D.M. 2018. Predicting habitat suitability for eleven imperiled fluvial freshwater mussels. Hydrobiologia 809(1): 265–283. DFO. 2007a. Documenting Habitat Use of Species at Risk and Quantifying Habitat Quality. DFO Can. Sci. Advis. Sec. Sci. Advis. Rep. 2007/038. DFO. 2007b. Revised Protocol for Conducting Recovery Potential Assessements. DFO Can. Sci. Advis. Sec. Sci. Advis. Rep. 2007/039 Haag, W.R. 2012. North American freshwtaer mussels: natural history, ecology, and conservation. Cambridge University Press, New York, NY. 505 p. Haag, W.R. 2019. Reassessing enigmatic mussel feclines in the United States. Freshw. Mollusk Biol. Conserv. 22(2): 43–60. Haag, W.R., and Rypel, A.L. 2011. Growth and longevity in freshwater mussels: evolutionary and conservation implications. Biol. Rev. 86(1): 225–247. Haag, W.R., and Staton, J.L. 2003. Variation in fecundity and other reproductive traits in freshwater mussels. Freshw. Biol. 48(12): 2118–2130. Haggerty, T.M., Garner, J.T., Patterson, G.H., and Jones, L.C. 1995. A quantitative assessment of the reproductive biology of Cyclonaias tuberculata (Bivalvia: Unionidae). Can. J. Zool. 73(1): 83–88. Hove, M. 1997. Ictalurids serve as suitable hosts for Purple Wartyback. Triann. Unionid Rep. 11(4). Jirka, K.J., and Neves, R.J. 1992. Reproductive biology of four species of freshwater mussels (Molluscs: Unionidae) in the New River, Virginia and West Virginia. J. Freshw. Ecol. 7(1): 35–44. Lamothe, K.A., van der Lee, A.S., Drake, D.A.R., and Koops, M.A. 2021. The translocation trade-off for eastern sand darter (Ammocrypta pellucida): balancing harm to source populations with the goal of re-establishment. Can. J. Fish. Aquat. Sci. 78(9): 1321–1331. Lande, R. 1988. Genetics and demography in biological conservation. Science 241(4872): 1455–1460. van der Lee, A.S., and Koops, M.A. 2021. Recovery Potential Modelling of Pygmy Whitefish (Prosopium coulterii) in Canada (Great Lakes – Upper St . Lawrence populations). Can. Sci. Advis. Sec. Res. Doc. 2021/026: iv + 20. McGowan, C.P., Runge, M.C., and Larson, M.A. 2011. Incorporating parametric uncertainty into population viability analysis models. Biol. Conserv. 144(5): 1400–1408. Metcalfe-Smith, J.L., McGoldrick, D.J., Zanatta, D.T., and Grapentine, L.C. 2007. Development of a monitoring program for tracking the recovery of endangered freshwater mussels in the Sydenham River, Ontario. Water Science and Technology Directorate, Environment Canada, Burlington, ON. 07-210: 63 p Morris, W.F., and Doak, D.F. 2002. Quantitative conservation biology: theory and practice of population viability analysis. Sinauer Associates, Sunderland, MA. 480 p. Palmqvist, E., and Lundberg, P. 1998. Population extinctions in correlated environments. Oikos 83(2): 359–367. Patterson, C.G. 1985. Biomass and production of the Unionid, Elliptio complanata (Lightfoot) in an old reservoir in New Brunswick. Freshw. Invertebr. Biol. 4(4): 201–207. https://www.dfo-mpo.gc.ca/csas-sccs/Publications/SAR-AS/2007/2007_038-eng.html https://www.dfo-mpo.gc.ca/csas-sccs/Publications/SAR-AS/2007/2007_039-eng.html https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2021/2021_026-eng.html https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2021/2021_026-eng.html https://publications.gc.ca/site/eng/9.868192/publication.html https://publications.gc.ca/site/eng/9.868192/publication.html https://publications.gc.ca/site/eng/9.868192/publication.html 24 R Core Team. 2021. R: a language and environment for statistical computing. Vienna, Austria. Reed, D.H. 2004. Extinction risk in fragmented habitats. Anim. Conserv. 7(2): 181–191. Reed, D.H., O’Grady, J.J., Ballou, J.D., and Frankham, R. 2003. The frequency and severity of catastrophic die-offs in vertebrates. Anim. Conserv. 6(2): 109–114. Roberts, J.H., Angermeier, P.L., and Anderson, G.B. 2016. Population viability analysis for endangered Roanoke logperch. J. Fish Wildl. Manage. 7(1): 46–64. Shaffer, M.L. 1981. Minimum population sizes for species conservation. Bioscience 31(2): 131– 134. Sheldon, M.N., McNichols-O’Rourke, K.A., and Morris, T.J. 2020. Summary of initial surveys at index stations for long-term monitoring of freshwater mussels in southwestern Ontario between 2007 and 2018. Can. Manuscr. Rep. Fish. Aquat. Sci. 3203: vii + 85 p. Sietman, B.E., Davis, J.M., and Hove, M.C. 2012. Mantle display and glochidia release behaviors of five quadruline freshwater mussel species (Bivalvia: Unionidae). Am. Malacol. Bull. 30(1): 39–46. Smith, M.W., Then, A.Y., Wor, C., Ralph, G., Pollock, K.H., and Hoenig, J.M. 2012. Recommendations for catch-curve analysis. North Am. J. Fish. Manag. 32(5): 956–967. Vélez-Espino, L.A., and Koops, M.A. 2009. Quantifying allowable harm in species at risk: application to the Laurentian black redhorse (Moxostoma duquesnei). Aquat. Conserv. Mar. Freshw. Ecosyst. 19(6): 676–688. Vélez-Espino, L.A., and Koops, M.A. 2012. Capacity for increase, compensatory reserves, and catastrophes as determinants of minimum viable population in freshwater fishes. Ecol. Model. 247: 319–326. Vélez-Espino, L.A., Randall, R.G., and Koops, M.A. 2010. Quantifying habitat requirements of four freshwater species at risk in Canada: Northern Madtom, Spotted Gar, Lake Chubsucker, and Pugnose Shiner. Can. Sci. Advis. Secr. Res. Doc. 2009/115. iv + 21 p. Young, J.A.M., and Koops, M.A. 2013. Recovery potential modelling of Hickorynut (Obovaria olivaria) in Canada. DFO Can. Sci. Advis. Sec. Res. Doc. 2013/022. v + 13 p. https://www.r-project.org/ https://publications.gc.ca/site/eng/9.889741/publication.html https://publications.gc.ca/site/eng/9.889741/publication.html https://publications.gc.ca/site/eng/9.889741/publication.html https://www.dfo-mpo.gc.ca/csas-sccs/publications/resdocs-docrech/2009/2009_115-eng.htm https://www.dfo-mpo.gc.ca/csas-sccs/publications/resdocs-docrech/2009/2009_115-eng.htm https://www.dfo-mpo.gc.ca/csas-sccs/publications/resdocs-docrech/2009/2009_115-eng.htm https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2013/2013_022-eng.html https://www.dfo-mpo.gc.ca/csas-sccs/Publications/ResDocs-DocRech/2013/2013_022-eng.html ABSTRACT INTRODUCTION METHODS SOURCES THE POPULATION MODEL PARAMETERIZATION Life-History Density-Dependence Stochasticity IMPACT OF HARM Elasticity of λ Simulation RECOVERY TARGETS Abundance: Minimum Viable Population (MVP) Habitat: Minimum Area for Population Viability (MAPV) RECOVERY TIMES RESULTS IMPACT OF HARM Elasticity of λ Simulation RECOVERY TARGETS Abundance: Minimum Viable Population (MVP) Habitat: Minimum Area for Population Viability (MAPV) RECOVERY TIMES DISCUSSION UNCERTAINTIES ELEMENTS REFERENCES CITED