\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Hunting alters viral transmission and evolution}



\author[1]{Nicholas Fountain-Jones}%
\author[2]{Simona Kraberger}%
\author[2]{Roderick Gagne}%
\author[3]{Marie Gilbertson}%
\author[2]{Daryl Trumbo}%
\author[1]{Michael Charelston}%
\author[2]{Patricia Salerno}%
\author[2]{Chris Funk}%
\author[2]{Kevin Crooks}%
\author[4]{Ken Logan}%
\author[4]{Mathew Alldredge}%
\author[5]{Simon Dellicour}%
\author[5]{Guy Baele}%
\author[6]{Xavier Didelot}%
\author[2]{Sue VandeWoude}%
\author[1]{Scott Carver}%
\author[7]{Meggan Craft}%
\affil[1]{University of Tasmania}%
\affil[2]{Colorado State University}%
\affil[3]{University of Minnesota Twin Cities}%
\affil[4]{Colorado Parks and Wildlife}%
\affil[5]{KU Leuven}%
\affil[6]{University of Warwick}%
\affil[7]{University of Minnesota}%


\vspace{-1em}



  
  \date{May 27, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
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 (higher purifying and diversifying selection)
compared with when hunting was reinstated. Changes in transmission could
be linked to short term social changes as male population increases.
These findings 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 profound impacts on pathogen
transmission and evolution not previously considered.%
\end{abstract}%



\sloppy


\textbf{Hunting alters viral transmission and evolution}

\textbf{Authors:} Nicholas M. Fountain-Jones \textsuperscript{1, 2*},
Simona Kraberger (simona.kraberger@gmail.com) \textsuperscript{3},
Roderick Gagne (rgagne@vet.upenn.edu) \textsuperscript{3}, Marie L.J.
Gilbertson (mjones029@gmail.com)\textsuperscript{1}, Daryl R. Trumbo
(daryl.trumbo@gmail.com) \textsuperscript{4}, Michael Charleston
(Michael.Charleston@utas.edu.au) \textsuperscript{2}, Patricia Salerno
(patriciasalerno@gmail.com) \textsuperscript{4,5}, W. Chris Funk
(Chris.Funk@colostate.edu) \textsuperscript{4}, Kevin Crooks
(kevin.crooks@colostate.edu) \textsuperscript{6}, Kenneth Logan
(ken.logan@state.co.us)\textsuperscript{7}, Mathew Alldredge
(Mat.alldredge@state.co.us) \textsuperscript{8}, Simon Dellicour
(simon.dellicour@kuleuven.be) \textsuperscript{9,10}, Guy Baele
(guy.baele@kuleuven.be) \textsuperscript{10}, Xavier Didelot
(xavier.didelot@warwick.ac.uk) \textsuperscript{11}, Sue VandeWoude
(Sue.Vandewoude@colostate.edu ) \textsuperscript{3}, Scott Carver
(Scott.Carver@utas.edu.au) \textsuperscript{2}, and Meggan E. Craft
(craft004@umn.edu) \textsuperscript{1,12}

\textbf{Affiliations}

\textsuperscript{1} Department of Veterinary Population Medicine,
University of Minnesota, St Paul, Minnesota 55108.

\textsuperscript{2} School of Natural Sciences, University of Tasmania,
Hobart Australia 7001.

\textsuperscript{3} Department of Microbiology, Immunology, and
Pathology, Colorado State University, Fort Collins, CO 80523.

\textsuperscript{4} Department of Biology, Graduate Degree Program in
Ecology, Colorado State University, Fort Collins, CO 80523.

\textsuperscript{5} Universidad Regional Amaz\selectlanguage{ngerman}ónica IKIAM, Km 7 Vía
Muyuna, Napo, Ecuador

\textsuperscript{6}Department of Fish, Wildlife, and Conservation
Biology, Colorado State University, Fort Collins, CO 80523.

\textsuperscript{7} Colorado Parks and Wildlife, Montrose, CO 81401 USA.

\textsuperscript{8} Colorado Parks and Wildlife, Fort Collins, CO 80526
USA.

\textsuperscript{9} Spatial Epidemiology Lab (SpELL), Université Libre
de Bruxelles, CP160/12 50, av. FD Roosevelt, 1050 Bruxelles, Belgium.

\textsuperscript{10} Department of Microbiology, Immunology and
Transplantation, Rega Institute, KU Leuven, Herestraat 49, 3000 Leuven,
Belgium.

\textsuperscript{11} School of Life Sciences and Department of
Statistics, University of Warwick, Coventry, CV47AL, United Kingdom.

\textsuperscript{12} Department of Ecology, Evolution and Behavior,
University of Minnesota, St Paul, Minnesota, 55108

*Nick.FountainJones@utas.edu\emph{.au}

\textbf{Key words:} Feline immunodeficiency virus, Mountain lion,
Selection, Perturbation effect, Phylogenetic diversity, Puma, \emph{Puma
concolor}

\textbf{Authorship statement:} NFJ conducted the analysis and wrote the
initial draft of the paper to which all authors contributed. KL and MA
studied the puma populations in the field and provided the blood
samples. SK, DT, PS, RG and SV collected virus and host genetic data.
SD, GB, MC and XD contributed to the phylogenetic and transmission tree
analyses. MG contributed to the spatial analysis. MEC, SV, KC, CF and SC
conceived of the project.

\textbf{Data accessibility:}

DNA sequences---GenBank accession: MN563193 - MN563239. All other data
and code to perform the analysis will be available on github.

\textbf{Running title: Hunting alters viral evolution}

\textbf{Article type:} Letter

\textbf{Number of words (abstract):} 144

\textbf{Number of words (body):} 4380

\textbf{Number of references:} 50

\textbf{Number of figures:} 3

\textbf{Number of tables:} 0

\textbf{Number of text boxes:} 0



\textbf{Abstract:}

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 (\emph{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 (higher purifying and diversifying
selection) compared with when hunting was reinstated. Changes in
transmission could be linked to short term social changes as male
population increases. These findings 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 profound
impacts on pathogen transmission and evolution not previously
considered.

Human actions commonly alter wildlife populations. A classic example of
an alteration is hunting, which often has density and demographic
effects on a population (Milner-Gulland \emph{et al.} 2003;
Whitman\emph{et al.} 2004; Packer \emph{et al.} 2009; Treves 2009).
However, the consequences of these actions on pathogen transmission and
evolution are largely unknown, and the few available studies report
contradictory findings. Theory predicts that for density-dependent
pathogens hunting-induced changes to density should decrease
transmission rates yet make little difference to transmission dynamics
for frequency-dependent pathogens. However, in practice models suggest
that reducing host density can decrease (Lloyd-Smith \emph{et al.} 2005;
Potapov \emph{et al.} 2012) or even increase pathogen transmission and
prevalence (Choisy \& Rohani 2006; Beeton \& McCallum 2011). 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 a `perturbation effect' (Woodroffe \emph{et al.} 2004,
2006; Carter \emph{et al.} 2007; Carr \emph{et al.} 2019). For example,
culling induced changes or `perturbations' to badger (\emph{Meles meles}
) territorial behavior was considered a driver of increased bovine
tuberculosis transmission between badgers and cattle (e.g.,
Woodroffe\emph{et al.} 2006). However, there is also evidence that
population reduction has little impact on canine rabies (Morters
\emph{et al.}2013) or Tasmanian devil facial tumor disease (Lachish
\emph{et al.}2010) dynamics. Recent advances in high-resolution pathogen
sequencing and analytic approaches can now elucidate patterns of
pathogen transmission and evolution (Smith \emph{et al.} 2015; Didelot
\emph{et al.} 2017; Grubaugh \emph{et al.} 2019) 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 (Logan \& Runge 2020b) 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 on pathogen-host
dynamics.

RNA viruses are ideal agents for examining the effect of hunting on
pathogen transmission and evolution. Genomic variation rapidly accrues
in RNA viruses, enabling estimation of fine-scale epidemiological
processes (such as transmission between hosts) and the basic
reproduction number R\textsubscript{0} (the average number of secondary
cases per infection, Biek \emph{et al.} 2015; Didelot \emph{et
al.}2017). 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
(Woolhouse\emph{et al.} 2014; Pybus \emph{et al.} 2015). For example, if
the cessation of hunting 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., Pybus \& Rambaut 2009). 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 (Fountain-Jones \emph{et al.} 2018). 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 (Webb 2000)). In contrast, phylogenetic clustering
(i.e., lower phylogenetic diversity than expected by chance (Webb
2000)), may be a marker of increased transmission events within a
population.

Here we leverage viral data collected from closely monitored puma
(\emph{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 immunodeficiency virus
(FIV\textsubscript{pco}), which is a host-specific pathogen considered
relatively benign and not associated with overt disease outcomes
(Biek\emph{et al.} 2003). Even though FIV\textsubscript{pco} is endemic
in puma populations, novel infections can spread in susceptible and
previously infected individuals (Malmberg \emph{et al.} 2019). Evidence
suggests FIV\textsubscript{pco} is transmitted via aggressive
interactions in most instances, although vertical transmission is also
possible (Biek \emph{et al.} 2003; Fountain-Jones \emph{et al.} 2017).
We analyzed these viral data in both regions using a transmission
network approach (Didelot \emph{et al.} 2017; Fountain-Jones \emph{et
al.} 2018) that incorporates a stochastic epidemiological model with
pathogen genomic data to trace transmission between individual puma. The
treatment region consisted of puma in a
\textasciitilde{}12000km\textsuperscript{2} area in western Colorado in
which hunting prior to our study was common practice (see Logan \& Runge
2020a). Hunting was excluded for a five-year period (2004-2009,
``no-hunting period'') and reinstated for a further five years
afterwards (2009-2014, ``hunting period''). During the no-hunting period
in the treatment region, the population of independent pumas (i.e.,
adults and sub-adults) increased from an estimated 23 (2005) to 57
(2009) individuals with much of this growth occurring 2007-2010 (Logan
\& Runge 2020a, i.e., after a two year lag 2004-6 hereafter ``Lag 1'').
Adult and sub-adult male survival was significantly higher in the
no-hunting period (Logan \& Runge 2020a). When hunting resumed, the
overall population declined after a lag of two years (hereafter `Lag
2'). However, the decline in abundance and demographic effects on males
were severe and rapid with males \textgreater{} 6 years old apparently
eliminated from the population after two hunting seasons (Logan \& Runge
2020b). In contrast, over the same 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. While
density was higher in the stable region (see Table S1), individual
survival probabilities in the stable region were unaltered across years
(Moss \emph{et al.} 2016). 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\textsubscript{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.

\textbf{Cessation of hunting shifts transmission networks and increases
R\textsubscript{0}}

We found that reducing hunting mortality had profound effects on
FIV\textsubscript{pco} transmission dynamics. Even though the
populations in the treatment and stable region were of comparable size
(Table S1), our estimates of R\textsubscript{0} for the same virus 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 significant, Fig. 1). Other model parameters,
such as generation time (time between initial
FIV\textsubscript{pco}infection and onward transmission, Fig. S2) and
the proportion of missing cases (Fig. S3) yielded similar estimates in
both regions. This burst of transmission in the treatment population was
likely a product of transmission between males as they were dominant in
the network (Fig. 1a). In the treatment population, males had an overall
mean weighted degree (i.e., the number of likely transmission events per
individual, weighted by probability of transmission occurring) double
that of females (0.14 compared to 0.37), only one putative transmission
event occurring between sexes, and no detected female-female
transmission events. When we assessed weighted degree homophily of
male-male transmission events (i.e., the number of edges only between
males), and simulations revealed that the dominance of male-male
transmission events in the network was not random (1000 simulated
annealing network iterations, \emph{p} \textless{} 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 (Logan \& Runge 2020b). Male survival rates
in the hunting period were also lower than for either sex in the
treatment region (Logan \& Runge 2020a). 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 between regions (0 vs 0.05), our simulations show that we could
not reject the null hypothesis that this difference was by chance
(\emph{p} = 0.692, Fig. S4b). Female-to-female transmission events
occurred between highly related females supporting our previous findings
of a significant role of host relatedness in FIV spread for puma living
in this region (Fountain-Jones \emph{et al.} 2019). Taken together, our
results indicate that lower hunting mortality was associated with an
increase in the number of transmission events and a transmission shift
towards males.

After hunting was prohibited, the greater survival and increasing
abundance of males probably resulted in greater competition between
males for mates. As the dominant transmission mode for
FIV\textsubscript{pco} is considered to be via aggressive contacts
(VandeWoude \& Apetrei 2006), increased male competition for mates
appears a probable explanation for the change 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 (\textasciitilde{}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 (Logan \& Sweanor 2001; Hornocker \& Negri 2010). 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 (Trumbo \emph{et
al.}2019). While our estimates suggest that we were able to sample
approximately 40\% of the FIV\textsubscript{pco} infections in both
regions (Fig. S3)---arguably good coverage for secretive, free-ranging
wildlife---our models account for this type of missing data
(Didelot\emph{et al.} 2017). For example, nearly all putative
transmission events we identified 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 difficult to observe.

\textbf{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 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 find
similar phylogenetic clustering in the stable region or hunting period
(Fig. 3a). Moreover, our transmission tree results and previous
phylogenetic reconstruction (Fountain-Jones \emph{et al.} 2019) 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 \emph{skygrowth} demographic
analyses (Volz \& Didelot 2018). Viral genetic diversity rapidly accrued
at the end of the no-hunting period (\textasciitilde{}2009) before
markedly declining after \textasciitilde{}2011 when hunting was
reinstated (Fig. 3b), closely mirroring male population size estimates
(\emph{R\textsuperscript{2}} = 0.8, \emph{p} = 0.010, Fig 3c). Female
population size was not significantly correlated to viral population
growth rate (\emph{R\textsuperscript{2}} = 0.190, \emph{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
have no data on puma behavior across time, it is possible that the
hunting-induced increase in male density increased competition for mates
and thus aggressive interactions (Logan \& Sweanor 2001). No such
increase in FIV\textsubscript{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 (\emph{p} = 0.01, likelihood ratio =
6.31). Purifying selection, potentially as a signature of rapid
transmission events (e.g., \emph{22} ), was dominant in both periods
(97.25\% sites \selectlanguage{greek}ω \selectlanguage{english}\textless{} 1), as is often the case in error-prone RNA
viruses, but stronger in the non-hunting period (\selectlanguage{greek}ω\selectlanguage{english}\textsubscript{2 nh} =
0, \selectlanguage{greek}ω\selectlanguage{english}\textsubscript{2h} = 0.1). In contrast, there was no shift in
evolutionary pressure in the same periods in the stable population
(\emph{p} = 0.5, likelihood ratio = 0.43). While impacting a smaller
proportion of the loci overall (2.79\% loci \selectlanguage{greek}ω \selectlanguage{english}\textgreater{} 1), there
was strong diversifying selection in the no-hunting period as well
(\selectlanguage{greek}ω\selectlanguage{english}\textsubscript{3 nh} = 21.46, \selectlanguage{greek}ω\selectlanguage{english}\textsubscript{3h} = 2.8). We
identified five FIV\textsubscript{pco} loci under diversifying selection
using the MEME routine in both regions (cutoff: \emph{p} [?] 0.1). Two of
these loci were only found in isolates in males and, based on our
transmission models, likely infected by FIV\textsubscript{pco} in the
no-hunting period. There was no signature of diversifying or purifying
selection in the envelope gene (\emph{env} ), which was surprising given
that \emph{env} is generally under greater evolutionary pressure as it
is responsible for the virus binding to the host cells (Kenyon \& Lever
2011). All loci under diversifying selection were detected in the
FIV\emph{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\textsubscript{pco} phenotype
occurring in this population. Systematic shifts in evolutionary pressure
are known to occur when viruses switch hosts (e.g., Tamuri \emph{et al.}
2009; Forni\emph{et al.} 2017); 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\textsubscript{pco} is largely apathogenic in puma and therefore our
findings demonstrate the types of changes in pathogen transmission
dynamics that can be caused by hunting induced changes in wildlife
populations.

\textbf{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. On the surface our results seem to contradict the perturbation
effect hypothesis (Carter \emph{et al.} 2007), as hunting in our
treatment region reduced the number of putative transmission events as
some theoretical models suggest (Lloyd-Smith \emph{et al.} 2005;
Potapov\emph{et al.} 2012). However, in our case the cessation of
hunting in a population (which was previously hunted) facilitated
demographic change via increased male survivorship and abundance that
perturbed the system to a different demographic state (Logan \& Runge
2020a). Even though the `perturbation' was different here, as reduced
hunting pressure may have resulted in more interactions, likely via
enhanced male-male competition, the underlying mechanism could be
similar. An expansion of the perturbation effect to include any
management action that leads 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\textsubscript{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 culling on disease dynamics in empirical systems
(Prentice \emph{et al.} 2019). 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
(Biek\emph{et al.} 2015). This molecular and analytical 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.

\textbf{Materials and Methods}

\emph{Study area and puma capture}

Our study was conducted in two regions in the Rocky Mountains in
Colorado separated by \textasciitilde{}500 km but at similar elevations
and with similar puma densities (Lewis \emph{et al.} 2015), vegetative
and landscape attributes, yet with differing degrees of urbanization
(see Fig. S10 and Lewis \emph{et al.} 2015). In the treatment region in
the Uncompahgre Plateau on the Western Slope of Colorado, blood samples
were taken from 114 individuals (2004-2011: 50 females \& 33 males,
2012-2014: 21 females \& 19 males) and monitored intensively (e.g., GPS
collars) until their death or the end of the study in 2014. In the
stable management region in the Front Range of Colorado, blood samples
were taken from 56 individuals from 2005-2014 (2005-2011: 11 females \&
5 males, 2012-2014: 21 females \& 19 males). Captured pumas were
anesthetized with established sedative and tranquilizer protocols (Logan
2012) and released after blood, serum, and oral swabs were collected.
Animal sex, age, and capture location were recorded. See Fountain-Jones
\emph{et al.} (2019) for sample storage, FIV\textsubscript{pco} DNA
extraction and sequencing details. In brief, for samples that were qPCR
positive for FIV\textsubscript{pco,} the complete \emph{ORFA} and
\emph{pol} gene regions were isolated using a nested PCR protocol
(Fountain-Jones \emph{et al.} 2019). Recombination was removed and the
genes were concatenated together. See Table S1 for a summary of the
sequence data and a comparison of study area size, host mortality, and
host genetic diversity between regions.

\emph{Transmission and phylogenetic trees}

We constructed transmission trees between pumas in each region using the
R package TransPhylo (Didelot \emph{et al.} 2017). TransPhylo uses a
time-stamped phylogeny to estimate a transmission tree to gain inference
into ``who infected whom'' and when. Briefly, 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) (Didelot \emph{et al.} 2017). This approach is
particularly useful for pathogens where the outbreak is ongoing, and not
all cases are sampled (Didelot \emph{et al.} 2017), as is the case here.
We leveraged our FIV\textsubscript{pco}Bayesian phylogenetic
reconstructions from previous work and focused on the two clades of
FIV\textsubscript{pco} that predominantly occurred in each region (see
Fountain-Jones \emph{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 (Didelot \emph{et al.} 2017). We assumed that
generation time could be drawn from a Gamma distribution (k = 2, \selectlanguage{greek}θ \selectlanguage{english}=
1.5) estimating onward transmission on average 3 years post-infection
(95\% interval: 0.3 - 8 years, based on average puma age estimates
(Logan \& Sweanor 2001)). Based on previous work (Lewis \emph{et al.}
2015), we were confident that the proportion of cases (\selectlanguage{greek}π\selectlanguage{english}) sampled was
high, therefore we set the starting estimate of \selectlanguage{greek}π \selectlanguage{english}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 \textgreater{} 200. We computed the posterior distributions of
R\textsubscript{0}, \selectlanguage{greek}π, \selectlanguage{english}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 \emph{igraph} package (Cs\selectlanguage{ngerman}árdi \&
Nepusz 2006) 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 \emph{igraph} .

\emph{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 fitted a
variety of probability distributions to edge weight and degree of both
treatment and stable regions, then used AIC to select the best fitting
target distribution. Edge density, network size and the number of
isolated nodes were fixed 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 male-male or female-to-female transmission
events) of the null networks to the observed and calculated a bootstrap
p-value.

\emph{Selection analyses}

To test if the demographic changes driven by hunting resulted in a
reduction in the intensity of natural selection on
FIV\textsubscript{pco} , we examined selective pressure in both time
periods in each region using the RELAX hypothesis testing framework
(Wertheim \emph{et al.} 2015). The method builds upon random effects
branch-site models (BS-REL) (Kosakovsky Pond\emph{et al.} 2011) that
estimates the \selectlanguage{greek}ω \selectlanguage{english}ratio (the ratio of non-synonymous to synonymous
mutations or dN/dS) along each branch from a discrete distribution of
three \selectlanguage{greek}ω \selectlanguage{english}ratio classes allowing selection pressure to vary across the
phylogeny (Wertheim \emph{et al.} 2015). A \selectlanguage{greek}ω \selectlanguage{english}ratio of one corresponds to
neutral selection with values \textgreater{} 1 being evidence for
diversifying (positive) selection along a branch, and \textless{} 1
evidence for purifying (negative) selection along a branch. Briefly,
RELAX tests for relaxation of selection pressure by dividing branches
into three subsets; test branches (T), reference branches (R) and
unclassified branches (U) (Wertheim \emph{et al.} 2015) with
\selectlanguage{greek}ω\selectlanguage{english}\emph{\textsubscript{T}} (resp. \selectlanguage{greek}ω\selectlanguage{english}\emph{\textsubscript{R}} ) being the
estimated dN/dS ratio on test (resp. reference) branches. The discrete
distribution of \selectlanguage{greek}ω \selectlanguage{english}is calculated using BS-REL for each branch class, and
then branches belonging to each subset are compared. The reference
estimates of \selectlanguage{greek}ω \selectlanguage{english}are raised to the power of \emph{k} (an intensity
parameter) so that ~\(\omega_{T}=\omega_{R}^{k}\) in order to simplify model
comparison. The null RELAX model is when the \selectlanguage{greek}ω \selectlanguage{english}distribution and thus
selective pressure is the same in R and T (when\emph{k} = 1). The null
model is compared to an alternate model (using a likelihood ratio test)
that allows \emph{k} to vary so that when\emph{k} \textgreater{}1
selection pressure on the test branches was intensified or \emph{k}
\textless{} 1 indicating that selection pressure has been relaxed
(Wertheim \emph{et al.} 2015). In the relaxed scenario, \emph{k}
\textless{} 1 branches in R are under stronger purifying and
diversifying selection compared to T branches (e.g., \selectlanguage{greek}ω \selectlanguage{english}shifts from 0.1
to 0.001 or from 10 to 2). See Wertheim \emph{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 classified as `U' as the majority had low phylogenetic
support (posterior probability \textless{} 0.6). To further interrogate
the sequence data to identify individual sites under selection, we
performed the MEME (mixed-effects model of evolution) pipeline (Murrell
\emph{et al.} 2012). We performed both MEME and RELAX models using the
Datamonkey web application (Weaver \emph{et al.} 2018).

\emph{Population growth rate}

We applied the non-parametric \emph{skygrowth} method (Volz \& Didelot
2018) to examine if the FIV\textsubscript{Pco} population growth rate
fluctuated 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 fitted
these models using MCMC (100,000 iterations) assuming that
FIV\textsubscript{Pco} population size fluctuated 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\textsubscript{Pco} population growth was related to male and female
population size estimates (Logan \& Runge 2020a). Measuring the
correlation between population size estimates and patterns of population
growth using generalized linear models (Gill \emph{et al.} 2016; Volz \&
Didelot 2018) was not feasible due to the relatively small size of this
dataset.

\emph{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, (Faith 1992)). 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 (Fountain-Jones \emph{et al.} 2019). To capture
phylogenetic uncertainty in these estimates, we utilized the
computational efficiency of the \emph{PhyloMeasures} R package algorithm
(Tsirogiannis \& Sandel 2016) to calculate SES.PD and apply this across
a 1000 tree subsample of posterior trees (Fountain-Jones \emph{et
al.}2019). Specifically, for each calculation of SES.PD we compared our
observed PD to a uniform null model (i.e., isolate samples are taken
with equal {[}uniform{]} probability). The code and data to perform
these operations as well as the transmission tree analysis above can be
found here: https://github.com/nfj1380/TransmissionTreeCode

\textbf{Acknowledgments:} This project was funded by the National
Science Foundation Ecology of Infectious Diseases research program
grants (DEB 1413925) and an Australian Research Council Discovery
Project Grant (DP190102020). M.L.J.G. was supported by the Office of the
Director, National Institutes of Health under award number NIH
T32OD010993. The content is solely the responsibility of the authors and
does not necessarily represent the official views of the National
Institutes of Health. SD is supported by the \emph{Fonds National de la
Recherche Scientifique} (FNRS, Belgium). GB acknowledges support from
the Interne Fondsen KU Leuven / Internal Funds KU Leuven under grant
agreement C14/18/094, and the Research Foundation -- Flanders (`Fonds
voor Wetenschappelijk Onderzoek -- Vlaanderen', G0E1420N). MEC was
funded by the National Science Foundation (DEB-1654609 and 2030509) and
the CVM Research Office UMN Ag Experiment Station General Ag Research
Funds. Any use of trade, firm, or product names is for descriptive
purposes only and does not imply endorsement by the U.S. Government.

\textbf{Acknowledgments:} This project was funded by the National
Science Foundation Ecology of Infectious Diseases research program
grants (DEB 1413925) and an Australian Research Council Discovery
Project Grant (DP190102020). M.L.J.G. was supported by the Office of the
Director, National Institutes of Health under award number NIH
T32OD010993. The content is solely the responsibility of the authors and
does not necessarily represent the official views of the National
Institutes of Health. SD is supported by the \emph{Fonds National de la
Recherche Scientifique} (FNRS, Belgium). GB acknowledges support from
the Interne Fondsen KU Leuven / Internal Funds KU Leuven under grant
agreement C14/18/094, and the Research Foundation -- Flanders (`Fonds
voor Wetenschappelijk Onderzoek -- Vlaanderen', G0E1420N). MEC was
funded by the National Science Foundation (DEB-1654609 and 2030509) and
the CVM Research Office UMN Ag Experiment Station General Ag Research
Funds. Any use of trade, firm, or product names is for descriptive
purposes only and does not imply endorsement by the U.S. Government.

\textbf{References:}

Beeton, N. \& McCallum, H. (2011). Models predict that culling is not a
feasible strategy to prevent extinction of Tasmanian devils from facial
tumour disease. \emph{J. Appl. Ecol.} , 48, 1315--1323.

Biek, R., Pybus, O.G., Lloyd-Smith, J.O. \& Didelot, X. (2015).
Measurably evolving pathogens in the genomic era. \emph{Trends Ecol.
Evol.}

Biek, R., Rodrigo, A.G., Holley, D., Drummond, A., Anderson, C.R., Ross,
H.A., \emph{et al.} (2003). Epidemiology, genetic diversity, and
evolution of endemic feline immunodeficiency virus in a population of
wild cougars. \emph{J. Virol.} , 77, 9578--9589.

Carr, A.N., Milleson, M.P., Hern\selectlanguage{ngerman}ández, F.A., Merrill, H.R., Avery, M.L.
\& Wisely, S.M. (2019). Wildlife management practices associated with
pathogen exposure in non-native wild pigs in Florida, U.S.\emph{Viruses}
, 11.

Carter, S.P., Delahay, R.J., Smith, G.C., Macdonald, D.W., Riordan, P.,
Etherington, T.R., \emph{et al.} (2007). Culling-induced social
perturbation in Eurasian badgers \emph{Meles meles} and the management
of TB in cattle: An analysis of a critical problem in applied
ecology.\emph{Proc. R. Soc. B Biol. Sci.} , 274, 2769--2777.

Choisy, M. \& Rohani, P. (2006). Harvesting can increase severity of
wildlife disease epidemics. \emph{Proc. R. Soc. B Biol. Sci.} , 273,
2025--2034.

Csárdi, G. \& Nepusz, T. (2006). The igraph software package for complex
network research. \emph{InterJournal Complex Syst.} , 1695.

Didelot, X., Fraser, C., Gardy, J. \& Colijn, C. (2017). Genomic
infectious disease epidemiology in partially sampled and ongoing
outbreaks. \emph{Mol. Biol. Evol.} , 34, msw075.

Faith, D.P. (1992). Conservation evaluation and phylogenetic
diversity.\emph{Biol. Conserv.} , 61, 1--10.

Forni, D., Cagliani, R., Clerici, M. \& Sironi, M. (2017). Molecular
Evolution of Human Coronavirus Genomes. \emph{Trends Microbiol.}

Fountain-Jones, N.M., Kraberger, S., Gagne, R., Trumbo, D.R., Salerno,
P., Funk, W.C., \emph{et al.} (2021). Host relatedness and landscape
connectivity shape pathogen spread in a large secretive
carnivore.\emph{Communications Biology} 4, 12
https://doi.org/10.1038/s42003-020-01548-2

Fountain-Jones, N.M., Packer, C., Troyer, J.L., VanderWaal, K.,
Robinson, S., Jacquot, M., \emph{et al.} (2017). Linking social and
spatial networks to viral community phylogenetics reveals
subtype-specific transmission dynamics in African lions. \emph{J. Anim.
Ecol.} , 86, 1469--1482.

Fountain-Jones, N.M., Pearse, W.D., Escobar, L.E., Alba-Casals, A.,
Carver, S., Davies, T.J., \emph{et al.} (2018). Towards an
eco-phylogenetic framework for infectious disease ecology. \emph{Biol.
Rev.} , 93, 950--970.

Gill, M.S., Lemey, P., Bennett, S.N., Biek, R. \& Suchard, M.A. (2016).
Understanding past population dynamics: Bayesian coalescent-based
modeling with covariates. \emph{Syst. Biol.} , 65, 1041--1056.

Grubaugh, N.D., Ladner, J.T., Lemey, P., Pybus, O.G., Rambaut, A.,
Holmes, E.C., \emph{et al.} (2019). Tracking virus outbreaks in the
twenty-first century. \emph{Nat. Microbiol.}

Hornocker, M.G. \& Negri, S. (2010). \selectlanguage{english}\emph{Cougar : ecology and
conservation} . University of Chicago Press.

Handcock M, Hunter D, Butts C, Goodreau S, Krivitsky P, Morris M (2018).
ergm: Fit, Simulate and Diagnose Exponential-Family Models for Networks.
The Statnet Project (http://www.statnet.org). R package version 3.9.4,
https://CRAN.R-project.org/package=ergm.

Karcher, M.D., Palacios, J.A., Bedford, T., Suchard, M.A. \& Minin, V.N.
(2016). Quantifying and mitigating the effect of preferential sampling
on phylodynamic inference. \emph{PLOS Comput. Biol.} , 12, e1004789.

Kenyon, J.C. \& Lever, A.M.L. (2011). The molecular biology of feline
immunodeficiency virus (FIV). \emph{Viruses} .

Kosakovsky Pond, S.L., Murrell, B., Fourment, M., Frost, S.D.W.,
Delport, W. \& Scheffler, K. (2011). A random effects branch-site model
for detecting episodic diversifying selection. \emph{Mol. Biol. Evol.} ,
28, 3033--3043.

Lachish, S., McCallum, H., Mann, D., Pukk, C.E. \& Jones, M.E. (2010).
Evaluation of selective culling of infected individuals to control
tasmanian devil facial tumor disease. \emph{Conserv. Biol.} , 24,
841--51.

Lewis, J.S., Logan, K.A., Alldredge, M.W., Bailey, L.L., VandeWoude, S.
\& Crooks, K.R. (2015). The effects of urbanization on population
density, occupancy, and detection probability of wild felids.\emph{Ecol.
Appl.} , 25, 1880--1895.

Lloyd-Smith, J.O., Cross, P.C., Briggs, C.J., Daugherty, M., Getz, W.M.,
Latto, J., \emph{et al.} (2005). Should we expect population thresholds
for wildlife disease? \emph{Trends Ecol. Evol.}

Logan, K. \& Runge, J. (2020a). Effects of hunting on a puma population
in Colorado. \emph{Wildl. Monogr.}

Logan, K.A. (2012). \emph{Wildlife Research Report, Colorado Division of
Parks and Wildlife:62} .

Logan, K.A. \& Runge, J.P. (2020b). \emph{Effects of hunting on a puma
population in Colorado} . Denver.

Logan, K.A. \& Sweanor, L.L. (2001). \selectlanguage{english}\emph{Desert puma : evolutionary
ecology and conservation of an enduring carnivore} . Island Press.

Malmberg, J.L., Lee, J.S., Gagne, R.B., Kraberger, S., Kechejian, S.,
Roelke, M., \emph{et al.} (2019). Altered lentiviral infection dynamics
follow genetic rescue of the Florida panther. \emph{Proc. R. Soc. B
Biol. Sci.} , 286, 20191689.

Milner-Gulland, E.J., Bukreeva, O.M., Coulson, T., Lushchekina, A.A.,
Kholodova, M. V., Bekenovil, A.B., \emph{et al.} (2003). Reproductive
collapse in saiga antelope harems. \emph{Nature} , 422, 135.

Morters, M.K., Restif, O., Hampson, K., Cleaveland, S., Wood, J.L.N. \&
Conlan, A.J.K. (2013). Evidence-based control of canine rabies: a
critical review of population density reduction. \emph{J. Anim. Ecol.} ,
82, 6--14.

Moss, W.E., Alldredge, M.W. \& Pauli, J.N. (2016). Quantifying risk and
resource use for a large carnivore in an expanding urban-wildland
interface. \emph{J. Appl. Ecol.} , 53, 371--378.

Murrell, B., Wertheim, J.O., Moola, S., Weighill, T., Scheffler, K. \&
Kosakovsky Pond, S.L. (2012). Detecting individual sites subject to
episodic diversifying selection. \emph{PLoS Genet.} , 8.

Packer, C., Kosmala, M., Cooley, H.S., Brink, H., Pintea, L., Garshelis,
D., \emph{et al.} (2009). Sport hunting, predator control and
conservation of large carnivores. \emph{PLoS One} , 4.

Potapov, A., Merrill, E. \& Lewis, M.A. (2012). Wildlife disease
elimination and density dependence. \emph{Proc. R. Soc. B Biol. Sci.} ,
279, 3139--3145.

Prentice, J.C., Fox, N.J., Hutchings, M.R., White, P.C.L., Davidson,
R.S. \& Marion, G. (2019). When to kill a cull: Factors affecting the
success of culling wildlife for disease control. \emph{J. R. Soc.
Interface} , 16.

Pybus, O.G. \& Rambaut, A. (2009). Evolutionary analysis of the dynamics
of viral infectious disease. \emph{Nat. Rev. Genet.}

Pybus, O.G., Tatem, A.J. \& Lemey, P. (2015). Virus evolution and
transmission in an ever more connected world. \emph{Proc. R. Soc. B
Biol. Sci.}

Smith, M.D., Wertheim, J.O., Weaver, S., Murrell, B., Scheffler, K. \&
Kosakovsky Pond, S.L. (2015). Less is more: An adaptive branch-site
random effects model for efficient detection of episodic diversifying
selection. \emph{Mol. Biol. Evol.} , 32, 1342--1353.

Tamuri, A.U., dos Reis, M., Hay, A.J. \& Goldstein, R.A. (2009).
Identifying Changes in Selective Constraints: Host Shifts in
Influenza.\emph{PLoS Comput. Biol.} , 5, e1000564.

Treves, A. (2009). Hunting for large carnivore conservation. \emph{J.
Appl. Ecol.}

Trumbo, D., Logan, K., Gagne, E., Kozakiewicz, C., Kraberger, S.,
Fountain-Jones, N., \emph{et al.} (2019). Urbanization impacts apex
predator gene flow but not genetic diversity across an urban-rural
divide. \emph{Mol. Ecol.} , 28, 4926--4940.

Tsirogiannis, C. \& Sandel, B. (2016). PhyloMeasures: a package for
computing phylogenetic biodiversity measures and their statistical
moments. \emph{Ecography (Cop.).} , 39, 709--714.

VandeWoude, S. \& Apetrei, C. (2006). Going wild: lessons from naturally
occurring T-lymphotropic lentiviruses. \emph{Clin. Microbiol. Rev.} ,
19, 728--762.

Volz, E.M. \& Didelot, X. (2018). Modeling the growth and decline of
pathogen effective population size provides insight into epidemic
dynamics and drivers of antimicrobial resistance. \emph{Syst. Biol.} ,
67, 719--728.

Weaver, S., Shank, S.D., Spielman, S.J., Li, M., Muse, S. V. \&
Kosakovsky Pond, S.L. (2018). Datamonkey 2.0: A modern web application
for characterizing selective and other evolutionary processes.\emph{Mol.
Biol. Evol.} , 35, 773--777.

Webb, C.O. (2000). Exploring the phylogenetic structure of ecological
communities: An example for rain forest trees. \emph{Am. Nat.} , 156,
145--155.

Wertheim, J.O., Murrell, B., Smith, M.D., Pond, S.L.K. \& Scheffler, K.
(2015). RELAX: Detecting relaxed selection in a phylogenetic
framework.\emph{Mol. Biol. Evol.} , 32, 820--832.

Whitman, K., Starfield, A.M., Quadling, H.S. \& Packer, C. (2004).
Sustainable trophy hunting of African lions. \emph{Nature} , 428,
175--178.

Woodroffe, R., Cleaveland, S., Courtenay, O., Laurenson, M.K. \& Artois,
M. (2004). Infectious disease in the management and conservation of wild
canids. In: \emph{The Biology and Conservation of Wild Canids} . Oxford
University Press, pp. 123--142.

Woodroffe, R., Donnelly, C.A., Jenkins, H.E., Johnston, W.T., Cox, D.R.,
Bourne, F.J., \emph{et al.} (2006). Culling and cattle controls
influence tuberculosis risk for badgers. \emph{Proc. Natl. Acad. Sci. U.
S. A.} , 103, 14713--14717.

Woolhouse, M.E.J., Adair, K. \& Brierley, L. (2014). RNA viruses: A case
study of the biology of emerging infectious diseases. In: \emph{One
Health} . American Society of Microbiology, pp. 83--97.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

\textbf{Fig. 1.} Males (blue nodes/puma silhouettes) were dominant in
the FIV\textsubscript{pco} transmission network in the treatment region
(a) whereas females (pink nodes/puma silhouettes) were more central in
the transmission network in the stable region (b). Nodes connected to
each other via edges indicate the probability of transmission in either
direction. Node size in the networks (left) is scaled based on the
number of edges estimated for each individual. Edge width is scaled
according to the probability of the transmission events, where wider
edges indicate a more likely transmission event (see Fig. S1).
R\textsubscript{0} estimates (with 95\% highest posterior density (HPD))
are based on the stochastic branching epidemiological model underlying
each transmission network (see \emph{Materials \& Methods} , (Didelot
\emph{et al.} 2017)). Transmission trees (right) show these putative
transmission events through time with branch color indicating how many
missing edges are likely between individuals. Orange: lag 1 period from
the start of the no-hunting period when males are recruited into the
study area and lag 2 period at the start of the hunting period until the
population declines; yellow: hunting pressure relieved; red: hunting
pressure resumed; grey: stable region. White nodes: unsampled
individuals estimate by the model. I: individuals that were likely
immigrants in this region based on (Trumbo \emph{et al.} 2019). See
Figure S2 for the FIV\textsubscript{pco} generation time distributions
for each region and Figure S3 for the estimate of missing cases across
years\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

\textbf{Fig. 2:} Infection time distributions from our transmission
network model for individuals involved in a putative transmission chain,
along with the likely direction of transmission (red arrows) and the
spatial context (see Fig. S5/S6 for information on other transmission
events in the treatment and stable region). Grey circles encompassing
puma silhouettes in the map insets represent known territorial overlap
between individuals (based on observations from K. Logan) and is not
representative of home range size. Orange: lag 1 period from the start
of the no-hunting period when males are recruited to the study area and
lag 2 period at the start of the hunting period until the population
declines, light yellow: hunting pressure relieved, red: hunting
pressure. Birth year is indicated by the cub silhouette, and death year
of M73 is indicated by the black horizontal line. The orange horizontal
line indicates when the FIV\textsubscript{pco} CO III lineage was
introduced into this population based on node estimates from
(Fountain-Jones \emph{et al.} 2019). Red horizontal lines indicate
transmission time distributions (overlap between infection time
distributions) and `trans' means `likely transmitted to'.

\par\null

\textbf{Fig. 3.} Eliminating hunting mortality led to: (a) overall
greater phylogenetic clustering (i.e., lower phylogenetic diversity) of
FIV\textsubscript{pco} isolates standardized for sample number and (b)
an increase in FIV\textsubscript{pco} population growth rate that was
(c) strongly correlated with male population size rather than (d) female
population size. (a) Standardized phylogenetic diversity (*: SES.PD,
standardized effect size phylogenetic diversity calculated from 1000
posterior trees) estimates revealed strong patterns of phylogenetic
clustering (smaller distances between isolates than expected by chance)
when hunting pressure was relieved (negative values of SES.PD).
Otherwise, FIV\textsubscript{pco} isolates were more dispersed across
the tree (SES. PD \textasciitilde{} 0, indicated by the dashed line).
Estimates of FIV\textsubscript{pco} prevalence (number of qPCR
positives/total number sampled) are provided next to each box and
whisker plot with number of individuals tested shown in parentheses (see
Fig. S8 for estimates of prevalence across years). (b) Viral population
growth rate was estimated using Bayesian phylodynamic reconstruction
(Volz \& Didelot 2018). See Fig. S7a for the corresponding skyline plot
(effective population size through time estimated via the\emph{phylodyn}
model (Karcher \emph{et al.} 2016)) for the treatment region and Fig.
S7b/c for complementary plots for the FIV\textsubscript{pco} clade
dominant in the stable region. Lg refers to the period between hunting
being reinstated and the start of population decline.

\selectlanguage{english}
\FloatBarrier
\end{document}

