Hunting alters viral transmission and evolution

Hunting can fundamentally alter wildlife population dynamics, but the consequences of hunting on pathogen transmission and evolution remain poorly understood. Here we present a study that leverages a unique landscape-scale experiment coupled with pathogen transmission tracing, network simulation and phylodynamics to provide insights into how hunting shapes viral dynamics in puma (Puma concolor). We show that removing hunting pressure enhances the role of males in transmission, increases the viral population growth rate and the role of evolutionary forces on the pathogen compared to when hunting was reinstated. Changes in transmission could be linked to short term social changes while the male population increased. These ndings are supported through comparison with a region with stable hunting management over the same time period. This study shows that routine wildlife management can have impacts on pathogen transmission and evolution not previously considered.


Main
Human actions commonly alter wildlife populations. A classic example of an alteration is hunting, which often has density and demographic effects on a population [1][2][3][4] . However, the consequences of these actions on pathogen transmission and evolution are largely unknown, and the few available studies report contradictory ndings. Theory predicts that for pathogens with density-dependent transmission, hunting-induced decreases in density should decrease transmission rates yet make little difference to transmission dynamics for frequency-dependent pathogens. In practice, empirical data and models suggest that reducing host density can either decrease 5,6 or even increase pathogen transmission and prevalence 7,8 . The complex interplay between host density, demography, and behavior also makes predicting the impacts of hunting on pathogen dynamics complex. Limited empirical work shows that population reduction can increase pathogen prevalence via social perturbation [9][10][11][12] . For example, cullinginduced changes or 'perturbations' to badger (Meles meles) territorial behavior was considered a driver of increased bovine tuberculosis transmission among badgers e.g., 9 . However, there is also evidence that population reduction has little impact on canine rabies 13 or Tasmanian devil facial tumor disease 14 dynamics. Recent advances in high-resolution pathogen sequencing and analytic approaches can now elucidate patterns of pathogen transmission and evolution [15][16][17] that were previously out of reach. Here we address the effects of hunting on pathogen dynamics by capitalizing on pathogen sequences collected from a detailed study on the demographic effects of hunting 18 as well as from sequences obtained over the same time period in a region where little hunting occurred. Our approach enables us to provide insights into the cascading consequences of hunting, and the cessation of hunting, on hostpathogen dynamics.
RNA viruses are ideal agents for examining the effect of hunting and the cessation of hunting on pathogen transmission and evolution. Genomic variation rapidly accrues in RNA viruses, enabling estimation of ne-scale epidemiological processes (such as transmission between hosts) and the basic reproduction number (R 0 ) 16,19 (see Box 1 for de nitions of key terminology highlighted in bold). Altered transmission dynamics and the arrival of new lineages can imprint distinctive evolutionary signatures on RNA viruses as they adapt quickly to changes in host populations they encounter 20,21 . For example, if a change of management led to a higher frequency of transmission events, we expect that the transmission bottleneck would lead to high purifying selection since within-host mutations are lost with transmission (e.g., 22 ). Conversely, if new mutations entering the host population allow the pathogen to escape immune detection, we may expect an increase in diversifying selection. Altered transmission dynamics and new lineages will also shape the phylogenetic diversity of the pathogen 23 . For example, if novel pathogen lineages are frequently arriving into a host population with limited transmission, we would expect to see a pattern of phylogenetic dispersion (i.e., higher phylogenetic diversity than expected by chance 24 ). In contrast, phylogenetic clustering (i.e., lower phylogenetic diversity than expected by chance 24 ) may be a marker of increased transmission events within a population.

Term
Definition R 0 The basic reproduction number 'R naught' is the expected number of cases generated by one case in a population of susceptible individuals. Transmission bottleneck Transmission of viruses between hosts usually involves a relatively small number of virus particles being exchanged between hosts (e.g., 53 ). This has the effect of reducing viral genetic diversity population size and creating a 'bottleneck'. Purifying selection 'Negative selection' is the removal of nonsynonymous mutations (i.e., mutations that lead to a change in protein coding). Diversifying selection 'Positive selection' is the favoring of nonsynonymous mutations that yield an adaptive advantage. These mutations can rapidly increase in frequency across a population. Transmission network A network where nodes represent individual puma and edges reflect transmission events based on transmission tree estimates. Edge weights are the probability of the transmission event occurring. Transmission trees generated by the R package 'TransPhylo' (Didelot et al, 2017) estimate who infected whom, including potentially unsampled individuals using a stochastic branching epidemiological model and a time-scaled phylogeny. Weighted degree The summed probability of a individual puma (i.e., a node in the network) being involved in transmission events divided by the number of transmission events (i.e., edges in the network). Weighted degree homophily The weighted degree of transmission events between members of the same sex.

Skygrowth demographic analyses
Non-parametric population-genetic model estimating the growth effective population size through time (a surrogate for genetic diversity) using Bayesian inference. This method has been shown to accurately reconstruct pathogen outbreak dynamics in a variety of systems ( 35,51 ).
Here we leverage viral data collected from closely monitored puma (Puma concolor) in two areas in Colorado during the same time period: a 'treatment region' in which hunting pressure changed over time and a 'stable management region' acting as a control (hereafter 'stable region'). We sequenced viral genes sampled from captured puma for an endemic RNA retrovirus, puma feline immunode ciency virus (FIV pco ), which is a host-speci c pathogen considered relatively benign and not associated with overt disease outcomes 25 . Even though FIV pco is endemic in puma populations, novel infections can spread in susceptible and previously infected individuals 26 . Evidence suggests FIV pco is often transmitted via aggressive interactions, although vertical transmission is also possible 25,27 . We analyzed these viral data in both regions using a transmission network approach 16,23 that incorporates a stochastic epidemiological model with pathogen genomic data to trace transmission between individual puma. The treatment region consisted of puma in a ~12000km 2 area in western Colorado in which hunting prior to our study was common practice (see 28  hereafter 'Lag 2'). However, the decline in abundance of males was severe and rapid with males > 6 years old apparently eliminated from the population after two hunting seasons 18 . In contrast, over the same 10-year period, the stable region in the Front Range of Colorado experienced minimal hunting pressures and no change in management practice. Nearly all the individuals sampled in both regions were adults and both sexes were evenly represented. Individual survival probabilities in the stable region were unaltered across years 29 . By comparing the treatment and stable regions, we were able to test how demographic changes caused by hunting cessation and reinstatement perturb viral transmission networks and epidemiological parameters (e.g., R 0 ), and also alter pathogen diversity and evolution. In doing so we begin to untangle the complex interplay between wildlife management and pathogen transmission, which is crucial for pathogen-orientated conservation and disease management strategies.

Cessation of hunting shifts transmission networks and increases R 0
We found that reducing hunting mortality had major effects on FIV pco transmission dynamics. Even though the populations in the treatment and stable region were of comparable size ( 41 , Table S1), our estimates of R 0 for the same virus over the 10-year period were two-fold higher in the treatment region compared to the stable region (with non-overlapping 95% high probability density intervals indicating that the difference is signi cant, Fig. 1). Other model parameters, such as generation time (the time between initial FIV pco infection and onward transmission, Fig. S2), and the proportion of missing cases (Fig. S3), yielded similar estimates in both regions. The burst of transmission in the treatment population after the cessation of hunting ( Fig. 1a right panel) was likely a result of transmission between males as they were dominant in the network. In the treatment population, males had an overall mean weighted degree (Box 1) double that of females (0.37 compared to 0.14). Only one putative transmission event occurred between sexes, and we detected no female-female transmission events. When we assessed weighted degree homophily of male-male transmission events, simulations revealed that the dominance of male-male transmission events in the network was not random (1000 simulated annealing network iterations, p < 0.001, Fig. S4a). Putative transmission events largely occurred when hunting mortality was eliminated ( Fig. 1a), during which time the survival of adults and subadult males was high, age structure increased, and the abundance of independent pumas increased 18 . Male survival rates in the hunting period were also lower than for either sex in the stable region 28 . Females were, however, much less connected in the transmission network in the treatment region compared to the stable region, where they were more central (Fig. 1b). In contrast to the treatment region, the stable region showed evidence of transmission from females to both females and males. Average weighted degree was higher overall for males than females in the stable region (0.46 vs 0.29). Even though weighted female-female degree homophily was higher in the stable region (0 vs 0.05), our simulations show that we could not reject the null hypothesis that this difference was by chance (p= 0.692, Fig. S4b). Female-to-female transmission events occurred between highly related females, supporting previous ndings of the importance of host relatedness in FIV pco spread for puma in this region 30 . Taken together, our results indicate that lower hunting mortality was associated with an increase in the number of transmission events which were dominated by males.
After hunting was prohibited, the greater survival and increasing abundance of males likely resulted in greater competition between males for mates. As the dominant transmission mode for FIV pco is considered to be via aggressive contacts 31 , increased male competition for mates appears a probable explanation for the differences in transmission dynamics. Further interrogation of our transmission network supports this theory, as in all but two instances, male-to-male transmission occurred between individuals with overlapping territories in the treatment region ( Fig. 2/S5/S6). One transmission pair was unusual in having less spatial proximity, yet one puma of this pair was a likely immigrant to the region (M133) and could have passed through M73's territory at some point (Fig. 2). With the exception of M73 (~6 y.o. at time of infection), all individuals involved in these transmission events were between 1-3 y.o., which is a period when males are establishing new territories and are starting to compete for access to females 32,33 . Our results suggest it is unlikely that these males transmitted to each other prior to dispersal or via maternal or paternal contacts-since these individuals were not related based on genomic data 34 . While our estimates suggest that we were able to sample approximately 40% of the FIV pco infections in both regions (Fig. S3)-arguably good coverage for secretive, free-ranging wildlifeour models account for this type of missing data 16 . For example, nearly all putative transmission events we identi ed from our transmission networks were between individuals on the landscape at the same time and in most cases were captured in close spatial proximity to each other. The biological plausibility of these transmission events demonstrates the power of adapting transmission network models to trace transmission and gain epidemiological insights in systems that are di cult to observe.
Hunting alters diversity and selective pressure on the virus Altered transmission dynamics at a population level were associated with changes in viral evolution and diversity in the treatment region. The increased number of transmission events in the no-hunting period compared to the hunting period was supported by the strong phylogenetic clustering (isolates with less phylogenetic diversity than expected by chance) detected relative to the hunting period (Fig. 3a). The link between reduced hunting pressure and increased transmission events was further supported as we did not nd similar phylogenetic clustering in the stable region or hunting period (Fig. 3a). Moreover, we found little evidence for new lineages arriving during the no-hunting period in the treatment region (Fig.  1a). We further interrogated viral diversity patterns across time using skygrowth demographic analyses 35 . Viral genetic diversity rapidly accrued at the end of the no-hunting period (~2009/2010) before markedly declining after ~2011 when hunting was reinstated (Fig. 3b), closely mirroring male population size estimates (R 2 = 0.8, p = 0.010, Fig 3c). Female population size was not signi cantly correlated to viral population growth rate (R 2 = 0.190, p = 0.630, Fig. 3d) adding further evidence for the enhanced role of male interactions in transmission dynamics when hunting mortality was reduced. While we lack behavioral observations of puma across time, it is possible that the increase in male density with the cessation of hunting allowed for increased competition for mates and thus aggressive interactions 33 . No such increase in FIV pco diversity and growth rate was detected in the stable population ( Fig.   S7b/c).
Within the treatment region, the increase in viral diversity was underpinned by greater effects of both purifying and diversifying selection acting on individuals infected during the no-hunting period compared to the hunting period (p = 0.01, likelihood ratio = 6.31). Purifying selection, potentially as a signature of rapid transmission events (e.g., 22 ), was dominant in both periods (97.25% sites ω < 1), as is often the case in error-prone RNA viruses, but stronger in the non-hunting period (ω 2 nh = 0, ω 2h = 0.1). In contrast, there was no shift in evolutionary pressure in the same periods in the stable population (p = 0.5, likelihood ratio = 0.43). While impacting a smaller proportion of the loci overall (2.79% loci ω > 1), there was strong diversifying selection in the no-hunting period as well (ω 3 nh = 21.46, ω 3h = 2.8). We identi ed ve FIV pco loci under diversifying selection using the MEME routine in both regions (cutoff: p ≤ 0.1). Two of these loci were only found in isolates in males and, based on our transmission models, the males were likely infected by FIV pco in the no-hunting period. There was no signature of diversifying or purifying selection in the envelope gene (env), which was surprising given that env is generally under greater evolutionary pressure as it is responsible for the virus binding to the host cells 36 . All loci under diversifying selection were detected in the FIV pol integrase region. Putting these lines of evidence together, we not only detected population-level impacts of demographic changes due to cessation of hunting on viral mutation, but also at the individual scale with stronger evolutionary pressure on viruses infecting males. Increased evolutionary pressure on the virus may increase the probability of a new FIV pco phenotype occurring in this population. Systematic shifts in evolutionary pressure are known to occur when viruses switch hosts e.g., 37,38 ; however here we show that selective constraints on a virus can be altered in response to host demographic changes caused by wildlife hunting. We stress that FIV pco is largely apathogenic in puma and therefore our ndings demonstrate the types of changes in pathogen transmission dynamics that can be caused by hunting induced changes in wildlife populations.
Perturbation, management and disease Our work provides a valuable case study on how hunting can have unexpected consequences for pathogen transmission and evolution across scales. Our multidisciplinary approach was particularly valuable in helping deconstruct how shifts in population structure and behaviour imprint on pathogen dynamics and evolution. For example, previous work using landscape genetic models only detected weak or inconsistent sex effects shaping FIV spread 27,30,39 . Our transmission network and phylodynamic approach, in contrast, was able to clearly distinguish the role of males in putative transmission chains and in accruing genetic diversity even though the data requirements are similar (e.g., a time scaled phylogeny). The putative transmission events we detected, supported by locational data, provided important mechanistic details at an individual scale that enabled us tease out the links between management, behaviour and transmission that are di cult to detect otherwise. Moreover, the shift in transmission network was able to provide context to the differences in pathogen evolution we detected between the no-hunting and no-hunting periods.
Our results provide a case study of the complex interplay between host demography, density and behaviour in shaping pathogen dynamics. In our case the cessation of hunting in a population in 2004 facilitated demographic change via increased male survivorship and abundance 28 , with potential increases in male-to-male contact behavior. Even though the 'perturbation' was the cessation of hunting, the underlying mechanism could be similar. An expansion in the way we think about perturbations to include a cessation of a practice leading to demographic or behavioral change may be warranted.
Our results also reveal potential shortcomings of relying on population estimates of prevalence to understand the impact of wildlife management actions on pathogen transmission. In our case, population estimates of FIV pco prevalence across time alone could not detect shifts in transmission associated with hunting and were not sensitive to changes in population size (Figs. S8/S9). The lack of signal from prevalence data may be a contributing factor behind the variability of the effects of hunting on disease dynamics in empirical systems 40 . Prevalence data may be better able to detect shifts in population demography where the pathogen causes acute infections with shorter periods of immunity.
The collection of pathogen molecular data from well-sampled wildlife populations across time is a logistical challenge, yet with ever cheaper and more mobile sequencing platforms, the potential to use approaches similar to ours is increasing, even for slowly evolving pathogens such as bacteria 19 . Our multidisciplinary approach can not only provide novel insights into the broader consequences of wildlife management on disease dynamics but can also help understand evolutionary relationships between hosts and pathogens in free-ranging species more broadly.

Materials And Methods
Study area and puma capture Our study was conducted in two regions in the Rocky Mountains in Colorado separated by ~500 km but at similar elevations and with similar puma densities 41 , vegetative and landscape attributes, yet with differing degrees of urbanization (see Fig. S10 Table S1 for a summary of the sequence data and a comparison of study area size, host mortality, and host genetic diversity between regions.

Transmission and phylogenetic trees
We constructed transmission trees between pumas in each region using the R package TransPhylo 16 . TransPhylo uses a time-stamped phylogeny to estimate a transmission tree to gain inference into "who infected whom" and when. Brie y, this approach computes the probability of an observed transmission tree given a phylogeny using a stochastic branching process epidemiological model; the space of possible transmission trees is sampled using reversible jump Markov chain Monte Carlo (MCMC) 16 . This approach is particularly useful for pathogens where the outbreak is ongoing, and not all cases are sampled 16 , as is the case here. We leveraged our FIV pco Bayesian phylogenetic reconstructions from previous work and focused on the two clades of FIV pco that predominantly occurred in each region (see Fountain-Jones et al. 2019). Whilst the TransPhylo approach makes few assumptions, a generation time distribution (the time from primary infection to onward transmission) is required to calibrate the epidemiological model 16 . We assumed that generation time could be drawn from a Gamma distribution (k = 2, θ = 1.5) estimating onward transmission on average 3 years post-infection (95% interval: 0.3 -8 years, based on average puma age estimates 33 ). Based on previous work 41 , we were con dent that the proportion of cases (π) sampled was high, therefore we set the starting estimate of π to be 0.6 (60% of cases tested in each region), and allowed it to be estimated by the model. We ran multiple MCMC analyses of 400,000 iterations and assessed convergence by checking that the parameter effective sample size (ESS) was > 200. We computed the posterior distributions of R 0 , π, and the realized generation time from the MCMC output. We also estimated likely infection time distributions for each individual and compared these estimates to approximate puma birth dates to ensure that these infection time distributions were biologically plausible. We then computed a consensus transmission tree for each region to visualize the transmission probabilities between individuals through time. Lastly, we reformatted the tree into a network object (nodes as individual puma and edges representing transmission probabilities) and plotted it using the igraph package 43 and overlaid puma sex as a trait.
Overall weighted degree and weighted degree for each sex, including edges representing homophily (e.g., male-male) and heterophily (e.g., male-female), were also calculated using igraph.

Simulation modelling
To test for non-random patterns of weighted degree between each sex, we applied a simulated network annealing approach from the Ergm R package (Handcock et al, 2018). To generate each simulated network, we tted a variety of probability distributions to edge weight and degree of both treatment and stable regions, then used AIC to select the best tting target distribution. Edge density, network size and the number of isolated nodes were xed based on each observed network. We added sex to each simulated node attribute drawing from a Bernoulli distribution (probability= 0.5). Using these network characteristics, we generated 1000 'null' networks and compared the homophily weighted degree distribution of each sex (i.e., the average weighted degree for each individual based on putative malemale or female-to-female transmission events) of the null networks to the observed and calculated a bootstrap p-value.

Selection analyses
To test if the demographic changes driven by hunting resulted in a reduction in the intensity of natural selection on FIV pco , we examined selective pressure in both time periods in each region using the RELAX hypothesis testing framework 44 . The method builds upon random effects branch-site models (BS-REL) 45 that estimates the ω ratio (the ratio of non-synonymous to synonymous mutations or dN/dS) along each branch from a discrete distribution of three ω ratio classes allowing selection pressure to vary across the phylogeny 44 . A ω ratio of one corresponds to neutral selection with values > 1 being evidence for diversifying (positive) selection along a branch, and < 1 evidence for purifying (negative) selection along a branch. Brie y, RELAX tests for relaxation of selection pressure by dividing branches into three subsets; test branches (T), reference branches (R) and unclassi ed branches (U) 44 with ω T (resp. ω R ) being the estimated dN/dS ratio on test (resp. reference) branches. The discrete distribution of ω is calculated using BS-REL for each branch class, and then branches belonging to each subset are compared. The reference estimates of ω are raised to the power of k (an intensity parameter) so that in order to simplify model comparison. The null RELAX model is when the ω distribution and thus selective pressure is the same in R and T (when k = 1). The null model is compared to an alternate model (using a likelihood ratio test) that allows k to vary so that when k >1 selection pressure on the test branches was intensi ed or k < 1 indicating that selection pressure has been relaxed 44 . In the relaxed scenario, k < 1 branches in R are under stronger purifying and diversifying selection compared to T branches (e.g., ω shifts from 0.1 to 0.001 or from 10 to 2). See Wertheim et al. (2015) for model details. T and R were selected from leaf branches (all other branches were Unassigned, U); individuals sampled from 2005-2011 (to the end of the lag period) were assigned to the R set and those sampled from 2012-2014 were assigned to T set. All branches not directly connecting to the tips were classi ed as 'U' as the majority had low phylogenetic support (posterior probability < 0.6). To further interrogate the sequence data to identify individual sites under selection, we performed the MEME (mixed-effects model of evolution) pipeline 46 . We performed both MEME and RELAX models using the Datamonkey web application 47 .

Population growth rate
We applied the non-parametric skygrowth method 35 to examine if the FIV Pco population growth rate uctuated across time and if this was related to changes in male or female population size in the treatment region. We did not do the same for the stable region as similar estimates were not available.
We tted these models using MCMC (100,000 iterations) assuming that FIV Pco population size uctuated every 6 months over a 14-year period (the estimated time to most recent common ancestor of this clade, Fig. S7). Otherwise, the default settings were used. We then performed a Pearson correlation test to assess if the trend in FIV Pco population growth was related to male and female population size estimates 28 . Measuring the correlation between population size estimates and patterns of population growth using generalized linear models 35,48 was not feasible due to the relatively small size of this dataset.

Phylogenetic diversity
To quantify phylogenetic diversity in each time period in each region, we calculated the standardized effect size (SES) for Faith's phylogenetic richness that accounts for differing sample sizes (SES for Faith's PD, 49 ). Faith's PD (hereafter PD) is the sum of the branch lengths of the phylogenetic tree linking all isolates for each subset (in this case the two time periods). As the number of isolates in each contrast differed (stable region 2005-2011: 11 isolates, stable region 2012-13: 5 isolates, treatment region 2005-2011: 10 isolates, treatment region 2012-14: 5 isolates) we calculated the standardized effect size (SES) by comparing the PD we observed to a null model that accounts for number of tips (i.e., how much phylogenetic diversity would we see for a given number of isolates by chance). We denote the standardized PD as SES.PD from here on; this was calculated across a subset of posterior phylogenetic trees from our previous Bayesian phylogenetic analyses 30 . To capture phylogenetic uncertainty in these estimates, we utilized the computational e ciency of the PhyloMeasures R package algorithm 50 to calculate SES.PD and apply this across a 1000 tree subsample of posterior trees 30  DNA sequences-GenBank accession: MN563193 -MN563239. All other data and code to perform the analysis will be available on github.  Figure S2 for the FIVpco generation time distributions for each region and Figure  S3 for the estimate of missing cases across year

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.