Assessing strategies to minimize unintended fitness consequences of aquaculture on wild populations

Artificial propagation programs focused on production, such as commercial aquaculture or forestry, entail strong domestication selection. Spillover from such programs can cause unintended fitness and demographic consequences for wild conspecifics. The range of possible management practices to minimize such consequences vary in their control of genetic and demographic processes. Here, we use a model of coupled genetic and demographic dynamics to evaluate alternative management approaches to minimizing unintended consequences of aquaculture escapees. We find that, if strong natural selection occurs between escape and reproduction, an extremely maladapted (i.e., nonlocal-origin, highly domesticated) stock could have fitness consequences analogous to a weakly diverged cultured stock; otherwise, wild population fitness declines with increasing maladaptation in the cultured stock. Reducing escapees through low-level leakage is more effective than reducing an analogous number of escapees from large, rare pulses. This result arises because low-level leakage leads to the continual lowering of wild population fitness and subsequent increased proportional contribution of maladapted cultured escapees to the total population. Increased sterilization efficacy can cause rapid, nonlinear reductions in unintended fitness consequences. Finally, sensitivity to the stage of escape indicates a need for improved monitoring data on how the number of escapees varies across life cycle stages.

A number of metrics can describe the genetic and demographic effects of escapees from cultivated populations on wild populations. In the main text we focus primarily on population size, fitness, and recovery time after cessation of an aquaculture program. Here we demonstrate that these metrics capture a wide variety of additional possibilities, with all metrics illustrated in the default case in fig. A.1. In addition, we test the sensitivity of our metrics to the run time chosen.
Two possible metrics of demographic effects of aquaculture are the equilibrium population size (fig. A.1a) and the fraction of natural spawners of natural origin (N W,t /N W,t+1 in eq. 7; fig. A.1g). Both show the maximum negative effect of aquaculture at the same value for the captive mean (θ C ) relative to the wild optimum. However, these metrics have opposing outcomes in comparing the extremes of similar versus different captive populations relative to the wild (θ C close to 1 versus 0, respectively). For more similar captive populations, more captive individuals survive natural selection, benefiting the wild population size and reducing the fraction of natural spawners of natural origin. For more different captive populations, captive individuals negatively affect the wild population through density-dependent mortality but then are less likely to survive natural selection, thus reducing the wild population size but resulting in a larger fraction of natural spawners being of natural origin.
For genetic effects of aquaculture, fitness (fig. A.1b, black solid line) is integrated over the full breeding value distribution (eq. 6) and thus captures aquaculture effects on both the genetic mean (μ W = gφ(g)dg; fig. A.1c) and variance (Ḡ W = (g −μ W ) 2φ (g)dg; fig. A.1d). Note that, as genetic variance changes with the evolution of the full breeding value distribution, heritability (h 2 =Ḡ W /(Ḡ W + E)) changes as well ( For the effect of aquaculture on the recovery of an impacted or endangered species, in the main text we focus on how long a wild population takes to recover after cessation of an aquaculture program ( fig. A.1e). A potentially relevant recovery question in the context of wild endangered species is whether the continued presence of an aquaculture program, with spillover, affects the recovery timeline for the wild population. We measure this potential recovery effect by the time it takes a small wild population (initializing the wild population size at N W,0 = 25 individuals) to recover to a threshold population size (defined as half of what the wild population would be at equilibrium without aquaculture) given aquaculture escapees ( fig. A.1f). Both recovery time metrics follow the same qualitative relationship with θ C , where the values of θ C that maximize recovery time without aquaculture lead to a wild population that remains below the threshold ("recovered") population size given continued aquaculture. Finally, in investigating the three central metrics with different model run times, we find that run time affects the quantitative outcome but not the qualitative trends ( fig. A.2). We also confirmed that the same comparative, qualitative trends across model runs as in figs. 1-3 and 5-7 in the main text occurred with a run time of t f = 30 time steps (time to ∼ 50% equilibrium in the default model scenario; results not shown).

B Breeder's equation simplification of the model
In order to better understand the dynamics of the strength of selection during the pulsed versus constant spillover simulations, we simplify the model to the analogous breeder's equation version. While providing a tractable metric of dynamical selection efficiency, this approach assumes a normal breeding value distribution with constant genetic variance and therefore underestimates the potential for aquaculture inputs to reduce fitness through increases in genetic variance. In the breeder's equation, the response to selection R is the product of the heritability h 2 and the selection differential S: The response to selection is the change in mean phenotypez t at time t, i.e., R = ∆z t = z t −z t . The heritability is the genetic proportion of the total phenotypic variance σ 2 z , which, assuming no genotype-by-environment interaction, is the sum of the additive genetic variance σ 2 g and environmental variance σ 2 e , i.e., h 2 = σ 2 g /σ 2 z where σ 2 z = σ 2 g + σ 2 e . The selection differential is the product of the phenotypic variance and the selection gradient β(z given mean fitnessW (z t ), i.e., S = σ 2 z β(z t ). Therefore, we can re-express eq. B.1 asz The mean fitnessW (z t ) is the product of the fitness W (z) of each phenotype z and the phenotypic probability distribution φ t (z) integrated over all phenotypes, i.e.,W (z t ) = W (z)φ t (z)dz. Given a normally-distributed phenotype distribution with meanz t and variance σ 2 and a fitness function representative of stabilizing selection for optimal trait θ with variance in the selection curve σ 2 w , then the mean fitness isW Using this function, the selection gradient and mean phenotype iteration are To apply the above dynamics to a population with aquaculture input, we first take separately the post-selectionz t (eq. B.7) for each of the wild (with mean phenotypez t ) and aquaculture (with mean phenotypez A ) populations. We then calculate the average phenotype of the combined population, weighted by relative input. This relative input is the post-selection size of each population, where we assume the mean fitness in eq. B.5 provides the origin-dependent survival. Given a reproductive output ofR and population size at time t of N t , the wild contribution isW (z t )RN t . Given a total aquaculture population size of N A and proportion release at time t of p t , the aquaculture contribution isW (z A )p t N A . Then the new mean phenotype is When following population size, we apply a maximum limit K such that This formulation assumes a life cycle ordering of reproduction-escape-selection-density dependence.
In order to determine selection efficiency at each point in time, we use the average selection differential (S t = σ 2 z β(z t ) given the selection gradient in eq. B.6) weighted by the relative contribution of each of the wild and aquaculture populations: To construct the figure presented in Box 1, we numerically iterate this model using analogous parameter values to those in the full model simulations from table 1: θ = θ W , σ 2 w = V S , σ 2 g = V LE , σ 2 e = V E ,R = ν I R, K = (R − 1)/(αR),z A = θ C , N C = N C , and p C = 0.05 in the constant-spillover simulations or, in the pulsed-spillover simulations, p C = 1 if the year is a multiple of twenty (0.05 −1 ) and 0 otherwise.

C Additional sensitivity results
Here we present the sensitivity of equilibrium population size and recovery time to the same parameter ranges explored in the fitness plot in the main text ( fig. 6).  fig. 6 for additional details on model implementation. Note differences in y-axis ranges.  fig. 3, first column, with escape after density dependence. Black lines indicate the outcome given constant, low-level spillover while gray lines give the outcome with stochastically variable spillover. For the pulsed spillover, the solid lines indicate the median outcome, dashed lines the 25th and 75th percentiles, and dotted lines the 1st and 99th percentiles of the data over multiple runs and years.
Here we present simulations equivalent to figs. 3 (constant versus pulsed escapement), 7 (ratio of captive escapees to wild-origin fish), and 8 (assortative mating) with escape occurring after density-dependent mortality, as opposed to before (the default assumption in the main text). with escape after density dependence. Different lines indicate simulations with different values for the mean aquaculture genotype (θ C , relative to the wild optimum phenotype of θ W = 1). The x-axis value of captive-origin:wild-origin population sizes is measured at escape.  fig. 8 with escape after density dependence. We implement assortative mating with increasing mating likelihood for increasing phenotypic similarity, where the parameter a represents the phenotypic correlation of mating pairs.