Predicting the Evolutionary Consequences of Trophy Hunting on a Quantitative Trait

Some ecologists suggest that trophy hunting (e.g. harvesting males with a desirable trait above a certain size) can lead to rapid phenotypic change, which has led to an ongoing discussion about evolutionary consequences of trophy hunting. Claims of rapid evolution come from the statistical analyses of data, with no examination of whether these results are theoretically plausible. We constructed simple quantitative genetic models to explore how a range of hunting scenarios affects the evolution of a trophy such as horn length. We show that trophy hunting does lead to trophy evolution defined as change in the mean breeding value of the trait. However, the fastest rates of phenotypic change attributable to trophy hunting via evolution that are theoretically possible under standard assumptions of quantitative genetics are 1 to 2 orders of magnitude slower than the fastest rates reported from statistical analyses. Our work suggests a re-evaluation of the likely evolutionary consequences of trophy hunting would be appropriate when setting policy. Our work does not consider the ethical or ecological consequences of trophy hunting.

). This has spawned the field of eco-evolution (Schoener, 2011). There is compelling empirical evidence 55 of rapid, joint phenotypic and ecological change from a number of systems (e.g., Hairston et al., 2005;Ozgul 56 et al., 2010), but evidence of genetic change is much less widespread (Yoshida et al., 2003), partly because it is 57 harder to demonstrate. Quite frequently, phenotypic change is attributable to evolution without supporting 58 evidence of genetic change (Hendry, 2016), or without examining whether the rates of evolutionary change 59 reported are theoretically plausible . argue that the activity has rapid detrimental consequences on hunted populations. However, no papers 67 have yet examined whether the rates of change observed by Coltman et al. (2003) and Pigeon et al. (2016) 68 are plausible using the quantitative genetic theory that motivated their statistical analyses, even though 69 skepticism has been raised as to whether the phenotypic changes observed can be attributed to evolution 70 (e.g., Traill et al., 2014). 71 We developed novel, general theory to examine the likely evolutionary consequences of selective harvesting 72 on a single sex in a sexually reproducing species. We worked in the quantitative genetics framework because 73 the genetic architecture of trophy traits is rarely known (e.g., Kruuk et al., 2002). We start with a brief 74 summary of quantitative genetic theory that motivated our models, and which is widely used to examine 75 the evolution of phenotypic traits of unknown genetic architecture in free-living populations (Merilä et al.,76 2001). We then describe the models we used, along with the parameter values we selected. 77

78
We use the following notation. Expectations and variances of the distribution of N (x, t) are denoted E(x, t) 79 and V (x, t) respectively. A subscript, either of f or m, is used to identify distributions or moments of 80 distributions taken over only females or males respectively. If this subscript is absent, the distribution is 81 taken over both sexes. We use a superscript R to identify distributions, or moments of distributions, that 82 have been operated on by selection. Quantitative genetics assumes that an individual's phenotype Z consists of the sum of various components.

85
These components include a breeding value A and the environmental component of the phenotype E, with 86 contributions from epistasis and non-additive genetic effects also sometimes included in the sum ( If alleles at a locus have an additive effect on a phenotypic trait, each allele can be assigned a value that 90 describes the contribution of that allele (in any genotype at that locus) to the phenotype. For example, 91 consider a bi-allelic locus with 3 genotypes, aa, aA and AA. Allele a has a value of 1g and allele A a value 92 of 2g. The breeding value of each genotype to body mass will be: aa = 1 + 1 = 2, aA = 1 + 2 = 3, and 93 AA = 2+2 = 4. Breeding values can be summed across genotypes at different loci to generate breeding values 94 for multi-locus genotypes. Under the additivity assumption, the dynamics of breeding values is identical to 95 the dynamics of alleles; this is the not always the case when non-additive genetic processes like heterozygote 96 advantage and epistasis are operating (Falconer, 1975).

97
Many applications of quantitative genetics use the infinitesimal model (Fisher, 1930). This assumes that 98 an individual's breeding value for a phenotypic trait is made up from independent contributions from a large 99 (technically infinite) number of additive genotypes, each making a very small contribution to the phenotypic 100 trait. There is no interaction between alleles at a locus (dominance) or interactions between genotypes at 101 different loci (epistasis).

102
In additive genetic models used to predict evolutionary change, it is usually assumed that E is determined by developmental noise. An individual's environmental component can be considered as a random value drawn from a Gaussian distribution with a mean and a constant variance: norm(0, V (E, t)). A and E are consequently independent. Thus, The distribution of breeding values is also assumed to be Gaussian, and These assumptions mean that, on average, the breeding value can be inferred from the phenotype -the 103 phenotypic gambit. The aim of statistical quantitative genetics is to correct the phenotype for nuisance 104 variables so the phenotypic gambit assumption is appropriate for the corrected phenotype (Lynch and Walsh, 105 1998).

106
Next, quantitative genetic theory makes the assumption that the mean of A among parents is equal to that in offspring: e.g., E(A, t + 1) = E R (A, t). In 2-sex models this requires that the expected value of A in an offspring is the mid-point of the breeding value of its parents. Given this assumption, where ∆E(Z) is the difference in the mean of the phenotype between the offspring and parental generations, 107 S(Z) is the selection differential on Z and V (A,t) V (Z,t) the heritability (h 2 ) of a trait, and t represents generation 108 number. The selection differential describes the difference in the mean value of the character between those 109 individuals selected to reproduce and the entire population prior to selection (Price, 1970). Equation () is 110 the univariate breeders equation (Falconer, 1975).

111
If all assumptions of the univariate breeders equation are met, it will accurately predict evolution of a 112 trait assuming that the selection differential and the additive genetic and phenotypic variances have been  (in the univariate case) and the G matrix (in the multivariate case).

125
A limitation of the breeders equation is it is not dynamically sufficient -it should not be used to make 126 predictions across multiple generations, particularly when evolution is sufficiently strong that it alters genetic 127 variances and covariances (Lande and Arnold, 1983). To construct a dynamic model, it is either necessary 128 to make assumptions about the genetic variance (it is sometimes assumed to be constant: (Lande, 1982 We developed a 2-sex, dynamic, quantitative genetic model to explore how hunting on one sex influences 133 phenotypic evolution. We iterate the population forwards on a per-generation time step. 134 We assume that in the absence of hunting, the trophy is not under selection in either sex and is con- are E(A, t) and E(E, t) at t = 0. Σ is a variance-covariance matrix, .

  
Variances can be chosen to determine the heritability h 2 at time t = 0.

141
We assume that males and females have the same distribution of phenotypes and breeding values at birth, 142 and that the birth sex ratio is unity: .

143
Next we impose selection. There is no direct selection on females and the number of recruits they 144 produced is set to 2, the replacement rate, to ensure the female population remains the same size over time We impose selection on males by culling a proportion α of individuals that are above average size, This generates a distribution of fathers N R m (A, E, t) that is equal in size to the distribution of mothers . 159 We now have distributions of maternal and parental characters that are the same sizes and sufficient 160 for the female population to replace itself with some males reproducing with multiple mothers. We assume

171
Taken together this gives the following recursion, Arnold, 1983). 179 We assume 2 traits Z 1 and Z 2 . We predict 1 generation ahead, so we do not use t for time to 184 We now impose selection on the phenotype with the following fitness function W (Z, t) = β 0 +β 1 Z 1 +β 2 Z 2 . 185 We estimate selection differentials on the 2 phenotypic traits as T is the vector transpose and S is a vector containing the selection differentials s 1 and s 2 . We also calculate 187 the univariate fitness functions W (Z) = β * 0 + β * 1 Z 1 and W (Z) = β 0 + β 2 Z 2 using methods from instrumental 188 variable analyses (Coulson et al., 2017;Kendall, 2015). From these functions, we calculated the univariate 189 selection differentials s * 1 and s 2 . We calculate univariate heritabilities using the relevant additive genetic  (Falconer, 1975).

196
To explore the effects of hunting on the evolution of a trophy, we ran simulations with a range of initial  We also examined the consequences of injecting additional genetic variance into the population at each

221
Selective trophy hunting led to an evolutionary response in all of our simulations ( Fig. 1-3). In our initial rate of loss of phenotypic variation and decline in the heritability scaled with harvesting rate (Fig. 1(B,C)).

230
When all males above the mean trophy value were harvested, additive genetic variance was initially rapidly 231 eroded, before starting to decline more slowly. This change was reflected in the dynamics of the phenotypic 232 variance ( Fig. 1(B)). These rates of change in the variance affected the dynamics of the mean phenotype.

233
Although the initial rate of evolution correlated with harvesting pressure, over the course of 100 generations 234 evolution was fastest when 75% of above average males were harvested. None of our scenarios predicted 235 phenotypic change at the rate reported by Coltman et al. (2003). In our initial simulations it took between 236 40 and 100 generations before the mean phenotype evolved to a value that would be significantly different 237 from its initial value (regardless of sample size). Finally, altering the initial heritability by reducing the initial 238 additive genetic variance slowed the rate of evolutionary changed as expected. In contrast, as the additive 239 genetic variance and consequently heritability increased, so too did the rate of evolution (Fig. 1(D)).

240
In our second simulation, we increased the initial heritability by reducing the environmental variation.

241
This had a relatively small impact on the rates of evolution (Fig. 2(A)), although the reduction in the 242 phenotypic variance (Fig. 2(B) did reduce rates of evolution at the highest levels of off-take (Fig. 2(A).

243
Increasing the additive genetic variance, and consequently the phenotypic variance, also increased rates of 244 evolutionary change slightly ( Fig. 3(A,B)), although rates of evolution were still between 1 and 2 orders

247
In all simulations, setting the segregation variance to a constant value generated linear selection because 248 selection does not rapidly erode the additive genetic variance (Fig. S2). However, even when all males of 249 above average horn size are culled, the rate of evolution is still > 5 times slower than that reported by breeders equation (Fig. 4(B)-(D)).

11
Although covariances in Σ(E) and Σ(A) could affect rates of evolution, when selection differentials were 265 large, covariances could not generate stasis or lead to evolution in the opposite direction to that predicted by 266 selection (Fig. 4, blue lines). However, as selection got weaker, correlated characters could prevent selection, 267 and even lead to very small evolutionary change in the opposite direction to that predicted by evolution (Fig.   268 4, red lines). However, effect sizes were small and would be challenging to detect without large quantities of regime that we assume.

291
A second assumption we make is that the trait is not subject to selection before selective harvesting  When predictions from simple models like ours fail to match with observation, the existence of genetically 307 correlated unmeasured characters is often invoked as an explanation (Merilä et al., 2001). Changing the 308 degree of generic covariation between two characters can significantly alter selection differentials on both 309 characters (e.g., Fig. 4). However, this does not mean that the failure to measure a correlated character will 310 lead to incorrect estimates of a selection differential on a trait. In fact, the failure to measure a correlated covariances that act to reduce the strength of selection can lead to low rates of evolutionary change in the 320 opposite direction to selection, but the effect is small and could only be detectable in very large data sets. 321 We consequently conclude that if the phenotypic gambit is assumed and significant selection on a trait is 322 observed, then unmeasured correlated characters can act to slow, or increase, rates of evolution compared to 323 those predicted by the univariate breeders equation, but they cannot result in evolutionary change that is 324 13 greater than the univariate selection differentials, or lead to evolutionary stasis. We conclude that although 325 our models on the effect of hunting on a trophy are simple, they will not be too wide of the mark, particularly 326 for large initial heritability for a trophy.

327
Although our models are simple, they provide some novel insights. In particular, our strongest selection 328 regimes result in initial increased rates of evolution. However, they erode the additive genetic covariance 329 more quickly than less stringent hunting regimes, rapidly slowing the rate of evolution. Over longer periods,

335
In most of our simulations we assume that the directional selection we impose erodes the additive genetic 336 variance as is often assumed in quantitative genetics (Falconer, 1975). We do this by constraining the have contributed to the observed phenotypic trends (Falconer, 1975). Quantitative genetics theory and em-360 pirical methods exists to deal with each of these processes (Lynch and Walsh, 1998), but statistical methods 361 to estimate these processes either require large population sizes or additional data that may not be avail- genotype-phenotype map, then neither will be the association between body size and horn length and fitness.

367
Although our models reveal that very rapid evolution attributable to selective hunting is not a plausible 368 explanation for the observed phenotypic declines, our models are not parameterized for bighorn sheep.  Quantitative genetics theory is powerful, elegant, and based on irrefutable logic (Falconer, 1975;Lande 381 and Arnold, 1983). The statistical methods used to estimate evolutionary change are also extremely powerful 382 when assumptions that underpin the analyses are met (Lynch and Walsh, 1998 Our work suggests that highly selective trophy hunting will result in evolutionary change, but that it will 397 not be particularly rapid. Evolutionary change would be more rapid if both sexes were selectively targeted lines, the greater the disparity between predictions from the univariate and multivariate breeders equation. 565 We simulated that all phenotypic variation is attributable to genetic (co)variances (A), approximately half 566 of phenotypic variance is attributable to additive genetic variance (B), and the effect of a positive (C) and 567 negative (D) covariance in the environmental components of the phenotypes on rates of evolution. The 568 genetic and environmental (co)variance used in each simulation can be found in Table S1.

569
Summary to the electronic given standard assumptions of quantitative genetics.