Apathogenic proxies for transmission dynamics of a fatal virus.

Identifying drivers of transmission-especially of emerging pathogens-is a formidable challenge for proactive disease management efforts. While close social interactions can be associated with microbial sharing between individuals, and thereby imply dynamics important for transmission, such associations can be obscured by the influences of factors such as shared diets or environments. Directly-transmitted viral agents, specifically those that are rapidly evolving such as many RNA viruses, can allow for high-resolution inference of transmission, and therefore hold promise for elucidating not only which individuals transmit to each other, but also drivers of those transmission events. Here, we tested a novel approach in the Florida panther, which is affected by several directly-transmitted feline retroviruses. We first inferred the transmission network for an apathogenic, directly-transmitted retrovirus, feline immunodeficiency virus (FIV), and then used exponential random graph models to determine drivers structuring this network. We then evaluated the utility of these drivers in predicting transmission of the analogously transmitted, pathogenic agent, feline leukemia virus (FeLV), and compared FIV-based predictions of outbreak dynamics against empirical FeLV outbreak data. FIV transmission was primarily driven by panther age class and distances between panther home range centroids. FIV-based modeling predicted FeLV dynamics similarly to common modeling approaches, but with evidence that FIV-based predictions captured the spatial structuring of the observed FeLV outbreak. While FIV-based predictions of FeLV transmission performed only marginally better than standard approaches, our results highlight the value of proactively identifying drivers of transmission-even based on analogously-transmitted, apathogenic agents-in order to predict transmission of emerging infectious agents. The identification of underlying drivers of transmission, such as through our workflow here, therefore holds promise for improving predictions of pathogen transmission in novel host populations, and could provide new strategies for proactive pathogen management in human and animal systems.

Introduction

Infectious disease outbreaks can have profound impacts on conservation, food security, and global health and economics. Mathematical models have proven a vital tool for understanding transmission dynamics of pathogens (1), but struggle to predict the dynamics of novel or emerging agents (2). This is at least partially due to the challenges associated with characterizing contacts relevant to transmission processes. Common modeling approaches that assume all hosts interact and transmit infections to the same degree ignore key drivers of transmission. Such drivers can include specific transmission-relevant behaviors including grooming or fighting in animals (3), concurrent sexual partnerships in humans (4), or homophily (5), and result in flawed epidemic predictions (6, 7). Further, identifying drivers of transmission and consequent control strategies for any given pathogen is typically done reactively or retrospectively in an effort to stop or prevent further outbreaks or spatial spread [e.g., (8, 9)]. These constraints limit the ability to perform prospective disease management planning tailored to a given target population, increasing the risk of potentially catastrophic pathogen outbreaks, as observed in humans (10), domestic animals (11), and species of conservation concern [e.g., (12–14)].

A handful of studies have evaluated whether common infectious agents present in the healthy animal microbiome or virome can indicate contacts between individuals that may translate to interactions promoting pathogen transmission (15–22). Such an approach circumvents some of the uncertainties associated with more traditional approaches to contact detection (6). In these cases, genetic evidence from the transmissible agent itself is used to define between-individual interactions for which contact was sufficient for transmission to occur. Results of such studies show mixed success (15–18). For example, members of the same household (19, 20) or animals with close social interactions (21, 22) have been found to share microbiota, but disentangling social mechanisms of this sharing is complicated by shared diets, environments, and behaviors (23).

These studies have, however, revealed ideal characteristics of non-disease inducing infectious agents (hereafter, apathogenic agents) for use as markers of transmission-relevant interactions. Such apathogenic agents should have rapid mutation rates to facilitate discernment of transmission relationships between individuals over time (24, 25). Furthermore, these agents should be relatively common and well-sampled in a target population, have a well-characterized mode of transmission that is similar to the pathogen of interest, and feature high strain alpha-diversity (local diversity) and high strain turnover (25, 26). RNA viruses align well with these characteristics (27) such that apathogenic RNA viruses could act as “proxies” of specific modes of transmission (i.e., direct transmission) and indicate which drivers underlie transmission processes. Such drivers, including but not limited to host demographics, relatedness, specific behaviors, or space use, could subsequently allow prediction of transmission dynamics of pathogenic agents with the same mode of transmission (25).

Here, we develop a novel workflow for identifying drivers of transmission in a naturally occurring host-pathogen system, and test the relevance of these drivers in the transmission of an analogously transmitted pathogenic virus. Florida panthers (Puma concolor coryi) are an endangered subspecies of puma found only in southern Florida. We have documented that this population is infected by several feline retroviruses relevant to our study questions (28, 29). Feline immunodeficiency virus (FIVpco; hereafter, FIV) occurs in ~50% of the population and does not appear to cause significant clinical disease (28). FIV is transmitted by close contact (i.e., fighting and biting), generally has a rapid mutation rate [intra-individual evolution rate of 0.00129 substitutions/site/year; (30)], and, as a chronic retroviral infection, can be persistently detected after the time of infection. Panthers are infected with feline leukemia virus (FeLV), also a retrovirus, which caused a well documented, high mortality outbreak among panthers in 2002–2004 (29). FeLV infrequently spills over into panthers following exposure to infected domestic cats (31, 32). Once spillover occurs, FeLV is transmitted between panthers by close contact and results in one of three infection states: progressive, regressive, or abortive infection (29). Progressive cases are infectious and result in mortality; regressive infections are unlikely to be infectious—though this is unclear in panthers—and recover (29, 33, 34). Abortive cases clear infection and are not themselves infectious (33).

The objectives of this study were therefore: (1) to determine which drivers shape FIV transmission in Florida panthers, and (2) test if these drivers are consistent with and can predict transmission of analogously transmitted FeLV in panthers. Success of this approach in our panther system would encourage testing similar apathogenic agents in other host-pathogen systems, with potential to improve our understanding of drivers of individual-level heterogeneity in transmission and consequently our ability to predict transmission dynamics of novel agents in human and animal populations.

Materials and methods

Dataset assembly

We assembled an extensive dataset covering almost 40 years of Florida panther research and including panther sex and age class. A subset of the population is monitored using very high frequency (VHF) telemetry collars, with relocations determined via aircraft typically three times per week. Previous panther research has generated a microsatellite dataset for monitored panthers (35), and a dataset of 60 full FIV genomes [proviral DNA sequenced within a tiled amplicon framework in (36)]. In addition, to augment observations from the 2002–04 FeLV outbreak (29), we leveraged an FeLV database which documents FeLV status (positive and negative) for 31 sampled panthers from 2002–04 as determined by qPCR.

FIV transmission inference

To determine drivers of FIV transmission, we first generated a “who transmitted to whom” transmission network using 60 panther FIV genomes collected from 1988 to 2011 [note that the panther population is small, with the average minimum annual panther counts across this period being 62.3 panthers; (37)]. We used the program Phyloscanner (38) (see Figure 1 for workflow across all analyses), which assumes both within- and between-host evolution when inferring transmission relationships between sampled and even unsampled hosts (38). Phyloscanner operates in a two step process, first inferring within- and between-host phylogenies in windows along the FIV genome. Then, using the within-host viral diversity gleaned from deep sequencing, Phyloscanner functionally performs ancestral state reconstruction to infer transmission relationships between hosts, outputting transmission trees or networks. For Phyloscanner's step one, we used 150bp windows, allowing 25bp overlap between windows. To test sensitivity to this choice, we separately ran a full Phyloscanner analysis with 150bp windows, but without overlap between windows (Supplementary methods). The tiled amplicon PCR approach used to generate our FIV genomic data biases for detection of one known variant, such that we did not expect detectable superinfections. In the second step of Phyloscanner, we therefore set the parameter which penalizes within-host diversity (k) to 0. We used a patristic distance threshold of 0.05 and allowed missing and more complex transmission relationships. Because we had uneven read depth across FIV genomes, we downsampled to a maximum of 200 reads per host. The output of the full Phyloscanner analysis was a single transmission network (hereafter, main FIV network).

Figure 1

Conceptual workflow across all analysis steps. Processes are shown on the left in blue; specific outcomes are shown on the right in green; the final analysis outcome is in yellow at the bottom right. Solid lines show direct flows or outcomes. Dashed lines show processes acting on or in concert with prior outcomes: for example, exponential random graph modeling (ERGM) was performed using the FIV transmission network, and the combination of the two produced the ERGM coefficients outcome.

To test sensitivity of our subsequent inference to variations in Phyloscanner output (e.g., due to the effects of random read downsampling, Phyloscanner windows, or sequencing errors), we also generated two summary FIV networks, varying the degree of window overlap in the first step of Phyloscanner analysis and re-running the random read downsampling in the second step. With Phyloscanner step one set to 25bp overlap, we generated four additional FIV transmission networks, but kept only those edges that were found in at least two of these four networks. We repeated this process with Phyloscanner step one set to 0bp overlap, again keeping only those edges found in at least two of four resulting transmission networks.

Statistical analysis of FIV transmission networks

Phyloscanner transmission tree output suggests direction of transmission, but in our case, the direction was often uncertain (see Results). To avoid putting undue emphasis on an uncertain direction of transmission, we simplified the transmission tree output to undirected, unweighted (binary) networks and performed statistical analysis of these networks using exponential random graph models [ERGMs; (39)]. ERGMs model the edges in networks, with explanatory variables representing the potential structural drivers of the observed network (39). By including network structural variables, ERGMs account for the inherent non-independence of network data. As such, we modeled “transmission relationships” (i.e., being connected in the transmission network) as a function of network structural variables and transmission variables we a priori expected to influence direct transmission processes in panthers. We considered several structural variables: an intercept-like edges term (39); geometrically weighted edgewise shared partner distribution (gwesp; representation of network triangles); alternating k-stars (altkstar; representation of star structures); and 2-paths [2 step paths from i to k via j; (40)]. In addition, we considered a suite of transmission variables (see Supplementary methods for additional variable details): panther sex; age class (subadult or adult); pairwise genetic relatedness [panther microsatellite data from (35)]; position of panther home range centroid (95% minimum convex polygon) or capture location (hereafter, centroid) relative to the major I-75 freeway (locations could be north or south of this east-west freeway); distance from centroid to nearest urban area [in km; USA Urban Areas layer, ArcGIS; (41)]; pairwise geographic distance between centroids (log-transformed; Supplementary Figure S1); and pairwise home range overlap [utilization distribution overlap indices of 95% bivariate normal home range kernels; (42, 43)]. We fit ERGMs for the main FIV network and the two summary FIV networks to verify robustness of inference.

Because ERGMs are prone to degeneracy with increasing complexity, we followed Silk and Fisher (39) and first performed forward selection for network structural variables, followed by forward selection of dyad-independent variables, while controlling for network structure. Model selection was based on AIC and goodness of fit, and MCMC diagnostics were assessed for the final model (Supplementary methods). ERGMs were fit with the ergm package (44) in R [v3.6.3; (45)].

Panther population and transmission simulations

We lack FeLV isolates to repeat Phyloscanner/ERGM analysis and thereby directly compare drivers of FeLV transmission to those identified for FIV. Rather, to determine the relevance of FIV transmission drivers for understanding and predicting FeLV transmission, we simulated FeLV transmission among panthers through a network structured by drivers of FIV transmission. We note that this approach is most representative of prospective disease modeling where models aim to predict transmission of a novel or emerging pathogen, and where transmission parameters are highly uncertain and models cannot be fit directly to data.

We first simulated panther populations that were representative of the population during the 2002-04 FeLV outbreak. Here, network edges represented likely transmission pathways based on ERGM-identified drivers of FIV transmission (FIV-based model). Hereafter, a full simulation includes both simulation of the panther population with its likely transmission pathways (i.e., a new network) and simulation of FeLV transmission within that population. This strategy of simulating new populations for each transmission simulation avoided putting excess weight on a small number of simulated and therefore uncertain networks. Below, we describe the process for a single simulation, but these procedures were repeated for each full simulation.

We first based the simulated population size on the range of empirical estimates from 2002–2004 [80–120 individuals; Supplementary Table S1; (37)]. Additional characteristics of the simulated population included those identified as significant variables in the ERGM analysis: age category and pairwise geographic distances between panther home range centroids (see Results). We randomly assigned age categories to the simulated population based on the proportion of adults vs. subadults. Age proportions were based on age distributions in the western United States (46), which qualitatively align with the historically elevated mean age of the Florida panther population [historically, mean age was as high as >6.5 years, but was about 4.5 years during 2002-04; (47)]. Pairwise geographic distances for the simulated population were generated by randomly assigning simulated home range centroids based on the distribution of observed centroids on the landscape (Supplementary methods).

We then used ERGM coefficients to generate network edges among the simulated panther population using the ergm package in R (44). The FIV transmission network spanned 15 years of observations and represents a subset of the actual contact network, as it includes only those interactions that resulted in successful transmission (48), and not non-transmission edges. We therefore had a high degree of uncertainty regarding the appropriate network density for our simulations. To manage and explore this uncertainty, we constrained density (ratio of existing edges to all possible edges) in our network simulations across a range of parameter space (net_dens, Supplementary Table S1).

The next step in each full simulation was to model FeLV transmission through the network generated from FIV predictors of transmission. FeLV transmission was based on a stochastic chain binomial process on the simulated network, following a modified SIR compartmental model (Figure 2). Simulations were initiated with one randomly selected infectious individual and proceeded in weekly time steps. Transmission simulations lasted until no infectious individuals remained or until 2.5 years, whichever came first.

Figure 2

Diagram of flows of individuals between compartments in the transmission model. Virus icons indicate infectious states, with the regressive infection icon darkened to represent reduced or uncertain infectiousness of this class. Note: a vaccination process was also included in the transmission model, but is not shown for simplicity. With vaccination, susceptibles could be vaccinated, and vaccinated individuals subsequently infected as with susceptibles, but with an additional probability of (1-ve). See Supplementary Table S1 for definitions of parameters.

Transmission was dependent on the following (Figure 2; see Supplementary Table S1 for parameter definitions): (1) existence of an edge between two individuals, (2) the dyad in question involving a susceptible and infectious individual, and (3) a random binomial draw based on the probability of transmission given contact (β). In addition, Puma concolor generally have low expected weekly contact rates (49); we therefore included a weekly contact probability, represented as a random binomial draw for contact in a given week (ω).

Upon successful transmission, infectious individuals were randomly assigned to one of three outcomes of FeLV infection (29). Progressive infections (probability P) are infectious (β), develop clinical disease, and die due to infection (μ). Regressive infections (also probability P) recover from infection (K*μ, where K is a constant ≤ 1) and, having entered a state of viral latency, are not considered at risk of FeLV reinfection (29, 50). Using model assumptions derived from known patterns of FeLV infection in domestic cats, regressive individuals are not infectious (29), but given ongoing uncertainty, we included some transmission from regressives (C*β, where C is a constant ≤ 1). Abortive infections (probability 1-2P) are never infectious, clearing infection and joining the recovered class. While the duration of immunity in abortive cases has not been studied in panthers, because abortive cases clear infection through a strong immune response and develop anti-FeLV antibodies, reinfection with FeLV is considered extremely unlikely (50).

A vaccination process was included in simulations as panthers were vaccinated against FeLV during the historical FeLV outbreak starting in 2003. Vaccination occurred at a rate, τ, and applied to the whole population, as wildlife managers are unlikely to know if a panther is susceptible at the time of capture or darting. However, only susceptible individuals transitioned to the vaccinated class (i.e., vaccination failed in non-susceptibles). Because panthers were vaccinated in the empirical outbreak with a domestic cat vaccine with unknown efficacy in panthers, we allowed vaccinated individuals to become infected in transmission simulations by including a binomial probability for vaccine failure (1-vaccine efficacy, ve, Supplementary Table S1).

The panther population size remained roughly static through the course of the FeLV outbreak (37). We therefore elected not to include background mortality, but did include infection-induced mortality. To maintain a consistent population size, we therefore included a birth/recruitment process. Because FIV-based simulated networks drew edges based on population characteristics, we treated births as a “respawning” process, in which territories vacated due to mortality were reoccupied by a new susceptible at rate, ν. This approach allowed us to maintain the ERGM-based network structure and is biologically reasonable, as vacated panther territories are unlikely to remain unoccupied for long. All simulations were programmed in R [v3.6.3; (45)].

Comparison of simulation predictions to observed FeLV outbreak

To evaluate the performance of our FIV-based model in the context of more common approaches used in predicting transmission of novel or emerging pathogens, we also predicted FeLV transmission dynamics using three alternative models: random networks, home range overlap-based networks, and a well-mixed model. The random networks model used Erdős-Rényi random networks, matching network densities from the FIV-based model (Supplementary Table S1), but otherwise allowing edges to occur between any pairs of individuals. Overlap-based networks were generated using the degree distributions of panther home range overlap networks from 2002 to 2004 and simulated annealing with the R package statnet [(51); Supplementary methods]. For both random and overlap-based networks, FeLV transmission was simulated as in the FIV-based simulations. The well-mixed model was a stochastic, continuous time compartmental model (Gillespie algorithm), with rate functions aligning with the chain binomial FeLV transmission probabilities (see Supplementary methods).

Consistent with modeling constraints when predicting transmission of novel or emerging agents, we performed transmission simulations for all model types (FIV-based, overlap-based, random, and well-mixed) across a range of reasonable parameter space (Supplementary Table S1), using a Latin hypercube design (LHS) to generate 150 parameter sets that efficiently sampled parameter space (52, 53). For each parameter set and model type, we performed 50 simulations (30,000 total). In each simulation, we recorded the number of mortalities and the duration of outbreaks, which were each summarized (medians) across each parameter set. To determine the accuracy of FIV-based predictions and this model's performance relative to alternative models, for each model type, we determined if each parameter set's predicted median (1) mortalities, (2) duration of outbreaks, and (3) abortive cases were within a reasonable range based on the observed FeLV outbreak [5–20 mortalities, 78–117 week duration, at least 5 abortive infections; (29)]. If so, a parameter set was deemed “feasible” for that model type. Ranges were used to account for uncertainty in observations and population size in this cryptic species (Supplementary methods). To compare the frequency of feasible FeLV predictions between model types, we fit a binomial generalized linear mixed model (GLMM), assuming a logistic regression with “feasible” (vs “unfeasible”) as the outcome, model type as a predictor variable, and a random intercept for parameter set.

We tested for spatial clustering of cases in the observed FeLV outbreak by leveraging our database of qPCR-based FeLV status. We performed a local spatial clustering analysis of FeLV cases and controls using SaTScan [50% maximum, circular window; (54)]. A SaTScan analysis seeks to identify clusters of cases in which the observed cases within a particular cluster exceed random expectation; this analysis reports the observed/expected ratio and radius of any significant clusters. In addition, we performed a global cluster analysis with Cuzick and Edward's test (global cluster detection with case-control data) in the R package smacpod [1, 3, 5, 7, 9, and 11 nearest neighbors; 999 iterations; (55–57)]. To determine if simulated FeLV cases demonstrated spatial clustering consistent with the observed outbreak, we repeated SaTScan local cluster analysis and Cuzick and Edward's tests (at 3, 5, and 7 nearest neighbors) with FIV-based simulation results. Because we would not expect representative spatial clustering in unfeasible parameter space (e.g., if epidemics were too large or small for spatial clustering to emerge), here alone we focused on the feasible subset of FIV-based simulation results. To verify that detected clustering in FIV-based simulations was not simply based on our respawning protocol, we also performed both spatial analyses with feasible overlap-based simulation results as a “negative control.” Because the overlap-based model was not spatially explicit, we assigned the same geographic locations to nodes in the overlap-based networks from the corresponding FIV-based networks.

To determine if feasible outcomes were especially sensitive to certain transmission parameters, we performed post hoc random forest variable importance analyses for each of the four model types with “feasible” vs. “unfeasible” as the binary response variable [using the R package randomForest (58, 59); see Supplementary results].

Results

FIV transmission network analysis

In the main FIV network, Phyloscanner inferred 42 potential transmission relationships (edges) between 19 individuals (nodes; network density = 0.25; Supplementary Table S2), after removing 9 edges that were between individuals known not to be alive at the same time (Figure 3). The summary transmission network allowing scanning window overlap included 20 nodes with 43 edges (network density = 0.23), and the summary network without window overlap included 20 nodes with 35 edges (network density = 0.18; after 8 and 6 edges removed, respectively, due to dates known alive). Panther FIV genomes missing from the transmission networks were those for whom transmission relationships could not be inferred by Phyloscanner (see Discussion).

Figure 3

Phyloscanner-derived main FIV transmission network. Node shape indicates panther age class (square = subadult; circle = adult). Node color indicates panther sex (blue = male; red = female). Edge weight represents Phyloscanner tree support for each edge (thicker edge = increased support); for visualization purposes, edges are displayed as the inverse of the absolute value of the log of these support values. While pictured as a directed and weighted network, statistical analyses used binary, undirected networks.

ERGM results for the main FIV network identified triangle (gwesp) and star structures (altkstar) as key structural variables, and age category and log transformed pairwise geographic distance as key transmission variables (Table 1; Supplementary Table S3). Though altkstar was not statistically significant, inclusion of this variable contributed to improved AIC and goodness of fit outcomes. Adults were more likely to be involved in transmission events (but see discussion of sample size limitations) and inferred transmission events were more likely between individuals which were geographically closer to each other. The fitted model showed reasonable goodness of fit (Supplementary Figures S2, S3). ERGM results for the two summary FIV transmission networks were comparable to the main FIV transmission network (Supplementary Table S3). The key difference was that the summary network with no window overlap did not find log-transformed pairwise geographic distances to be a significant variable, though this fitted model showed evidence of degeneracy. To further confirm consistency of our Phyloscanner and ERGM-based inference, we performed a post hoc analysis with simulated random networks, finding our results were generally robust to variations in Phyloscanner output (Supplementary results; Supplementary Figure S4).

Table 1

Main FIV transmission network exponential random graph model results.

Variable	Estimate	SE	p-value
Edges (intercept)	−2.56	1.33	0.055
gwesp	0.98	0.26	<0.001
altkstar	−0.70	0.96	0.47
Age (Adult)	0.93	0.44	0.03
Log pairwise distance	−0.45	0.21	0.03

“gwesp” is geometrically weighted edgewise shared partner distribution (a representation of triangle structures) and “altkstar” is alternating k-stars (a representation of star structures). Age classes were subadult and adult, with subadults the reference level; pairwise distances were between home range centroids and log-transformed. Only those variables from the final model are shown. Estimates shown are untransformed; SE represents standard error; p-values <0.05 were considered statistically significant.

FeLV simulations

About 9% of parameter sets across all model types were classified as feasible (Supplementary Figures S6, S7). The FIV-based model had the highest odds of feasibility, though this difference did not achieve statistical significance (Table 2). SaTScan analysis of observed FeLV status found weak evidence of local spatial clustering (two clusters detected, but not statistically significant with p = 0.165 and 0.997, respectively; Supplementary Figure S5). Cuzick and Edward's tests found evidence of global clustering at 3, 5, and 7 nearest neighbor levels (test statistic T where k is number of nearest neighbors considered: T3 = 20, p = 0.049; T5 = 32, p = 0.028; T7 = 43, p = 0.023). Feasible parameter sets from both the FIV-based and overlap-based models produced some evidence of local and global spatial clustering of simulated FeLV cases (Figure 4, Supplementary Figure S8). However, the FIV-based model better captured the size and strength of predicted local clusters (SaTScan radius and observed/expected cases, respectively; Figure 4) and was moderately better at capturing global spatial patterns (Supplementary Figure S8).

Table 2

Fixed effects results from model-type performance GLMM.

Variable	Estimate	SE	p-value
Intercept	0.055	0.40	<0.001
FIV-based network model	1.55	0.42	0.30
Random network model	1.32	0.43	0.52
Overlap-based network model	1.21	0.44	0.66

Estimates provided are exponentiated; the well-mixed model was the reference group and none of the model-type results achieved statistical significance.

Figure 4

SaTScan cluster analysis for feasible FIV-based and overlap-based network simulations show stronger agreement for the FIV-based model, compared to the overlap-based model, between empirical observations (red horizontal lines) and model predictions for (A) FeLV cluster size and (B) Observed/Expected FeLV cases associated with the top detected cluster. The overlap-based model, with locations assigned based on matching FIV-based simulations, served as a “negative control” for comparison to the FIV-based model's spatial predictions. Shown are feasible simulation results in which at least one cluster was detected with p-values ≤ 0.1; further, if SaTScan identified more than one cluster, only the results from the most well supported (i.e., top cluster) are shown.

The post-hoc random forest analyses typically showed poor balanced accuracy and area under the curve (AUC) results. However, the parameter shaping transmission from regressively infected individuals (C), consistently showed support for weak to moderate transmission from regressives (i.e., C = 0.1 or 0.5; Supplementary Figure S11).

Discussion

In this study we develop a new approach whereby we leverage genomic and network approaches to identify drivers of transmission of a common apathogenic agent. Further, we demonstrate that these drivers are relevant to and capable of prospectively predicting dynamics of an uncommon and virulent pathogen. Our approach was distinctly different from simpler models we tested, as the apathogenic (FIV)-based approach focused on underlying drivers or mechanisms of transmission and could be used to prospectively identify management-relevant predictors of transmission and develop disease control plans prior to an outbreak of a virulent pathogen (FeLV). We found that FIV transmission in panthers is primarily driven by adults and proximity between home range centroids, and that our FIV-based drivers of transmission predicted FeLV transmission dynamics at least as well as simpler alternative models in a prospective modeling framework (comparable to predicting transmission dynamics of novel or emerging pathogens). While we do not propose that this apathogenic agent approach could accurately predict exactly when, where, and to whom transmission might occur, our results support the role of apathogenic agents as novel tools for prospectively determining sources of individual-level heterogeneity in transmission and consequently improving proactive disease management.

FIV-based transmission drivers are relevant for FeLV transmission dynamics

We found that our network model based on drivers of FIV transmission produced FeLV outbreak predictions consistent with the observed FeLV outbreak. The FIV-based approach performed at least as well as simpler models, per our GLMM analysis, with evidence that FIV better predicted the observed spatial dynamics for FeLV transmission. A key difference between the FIV-based approach and other spatially explicit methods is that FIV allowed us to determine the importance of spatial dynamics prospectively (via transmission tree and ERGM analyses) and then translate to predictions of FeLV transmission, rather than relying on retrospective FeLV spatial analyses. Furthermore, while more complex potential drivers of transmission (e.g., host relatedness or assortative mixing by age or sex) were not found to be important for FIV transmission in this host-pathogen system, these may yet be key for driving transmission in other systems. Simpler model types like random networks or metapopulation models may struggle to make transmission predictions that incorporate these drivers of transmission-relevant contact. The predictive capabilities we observed here using drivers of an apathogenic virus could thus open new opportunities to determine behavioral and ecological drivers of individual-level heterogeneity in the context of pathogen transmission, and even shape epidemic management strategies for pathogens such as FeLV.

Our network statistical analysis (ERGMs) determined that pairwise geographic distances and age category structure FIV transmission in the Florida panther. These findings were generally robust to variations in the transmission network and are well supported by panther and FIV biology, providing confidence in the functioning of our workflow for identifying drivers of transmission. For example, panthers are wide-ranging animals but maintain home ranges, and this appears to translate to increased transmission between individuals that are close geographically. This finding is supported by the tendency for FIV phylogenies to show distinct broad (60) and fine scale (61) geographic clustering in Puma concolor. Further, specifically among Florida panthers, spatial autocorrelation of FIV exposure status was previously found to approach statistical significance (62). The wide-ranging nature of puma appears to limit geographic clustering of many infectious agents (62), with FIV a notable exception to this pattern. In addition, because FIV is a persistent infection, we would expect cumulative risk of transmission to increase over an individual's lifetime and adults would consequently be involved in more transmission events. The low number of subadult individuals in our dataset, however, means that this finding must be interpreted with some caution.

With these ERGM results in mind, key components of the success of our FIV-based approach are likely that (1) FIV is a largely species-specific virus with transmission pathways closely matching intraspecific transmission of FeLV, and (2) both FIV and FeLV, perhaps unusually for infectious agents of puma, display spatial clustering of infection. Here, FIV fundamentally acted as a proxy for close, direct contact in panthers, and could consequently determine drivers of such contacts. If, for example, FIV also exhibited strong vertical or environmental transmission, we would no longer expect the predictive success for FeLV we observed here. This consideration highlights the importance of careful apathogenic agent selection when attempting to identify drivers of transmission relevant to novel or emerging pathogen transmission. For example, the mixed results when using commensal agents to identify close social relationships in other systems (15–18, 21, 22) highlights that some host-apathogenic agent combinations will work better than others for determining drivers of transmission. Within our study, Phyloscanner struggled to elucidate transmission relationships between many of our FIV genomes, likely due to unusually low genetic diversity among our FIV isolates, or our use of proviral DNA (which has lower diversity than circulating RNA) (36). While the drivers of transmission we identified are biologically reasonable, we may have lacked the power to identify more complex relationships (e.g., homophily) due to the low number of individuals in our transmission network.

We propose that apathogenic agent selection should carefully consider agent genetic diversity within a target population—not just expected diversity based on typical mutation rates (24–27), as in our case—and favor those agents with high diversity to facilitate transmission inference. We also propose that apathogenic agents should represent the timescales of transmission for the pathogen of interest. For example, FeLV spreads slowly through panthers (31, 32), such that the transmission relationships identified by FIV, a chronic infectious agent that spreads among panthers over the course of years (36), may be most representative across the longer timescales we evaluated here. In contrast, short, acute pathogen epidemics would likely best be represented by apathogenic agent transmission over shorter timescales. Similarly, the timescale of data collection should correspond to the apathogenic agent of interest to reduce the probability of missing individuals in the inferred transmission network. Our results reinforce that, perhaps most importantly, an apathogenic agent should have a well characterized mode of transmission that closely matches transmission of the pathogen of interest (26), as this was likely key to our success with FIV and FeLV. Future research could determine how divergent an apathogenic agent may be from a pathogen of interest while still predicting transmission dynamics.

Potential applications

Our FIV-based approach to identifying drivers of transmission required extensive field sampling, though this is not infeasible in wildlife species of conservation concern or many livestock systems (27, 63). With increasing availability of virome data and even field-based sequencing technology, our proposed approach for identifying drivers of transmission relevant to predicting novel or emerging pathogen dynamics may become more accessible with time. From a practical perspective, if the only concern for prospective analysis of a pathogen of concern is predicting epidemic size and duration, our results indicate that a simpler approach would perform just as well as our approach. However, by identifying underlying drivers of transmission, our method also identified important, management-relevant spatial dynamics of transmission for FIV—dynamics which are uncommon among other infectious agents of wide ranging panthers (62). This presents a particular advantage to studies focused on identifying drivers of transmission, even among apathogenic agents [e.g., (21)]. For example, the addition of interventions such as vaccination to simulation models such as our FIV-model used here can help determine conditions (i.e., parameter space) in which spatially-targeted vaccination may be most effective (64, 65).

Further, our approach for identifying drivers of transmission could be applied in an adaptive management framework (66–68), in which apathogenic agent-based transmission predictions provide (1) a priori expectations for novel or emerging agent transmission dynamics that can aid in proactively designing targeted intervention strategies, and (2) a platform for updating strategies as new information becomes available in the event of an outbreak. Indeed, we have used a similar approach to determine optimal FeLV management strategies in panthers, including exploring a broad range of parameter space to determine how uncertainties in transmission parameters affect expected outcomes (69). While our transmission tree and ERGM results with FIV point to the role of spatial proximity for transmission, our method could similarly identify sex- or rank-biased transmission, homophily, or other transmission drivers relevant to pathogen management [e.g., (70, 71)]. We propose that our approach for identifying transmission drivers is best suited for proactive pathogen management in species of conservation concern, populations of high economic value (e.g. production animals), populations with infrequent pathogen outbreaks that make targeted surveillance more difficult, or populations at high risk of spillover (72), all of which may most benefit from rapid, efficient epidemic responses.

Caveats and future directions

While few parameter sets in our simulations were classified as feasible, this appears to be predominantly the result of the wide range of parameter space explored through our LHS sampling design. This limitation was fundamentally due to uncertainties in FeLV transmission parameters, and is representative of typical uncertainties experienced in predicting transmission of emerging or understudied pathogens (2, 73). Our method could also be applied in cases where the concern is introduction of a known pathogen into a new population (e.g., foreign animal diseases of livestock). In such cases, the parameter space to be explored may be much reduced [e.g., (74)]. Regardless of the confidence in transmission parameters, sensitivity analyses with variable importance analysis can highlight key parameters important for model outcomes [e.g., as in White et al. (59)]. If factored into an adaptive management plan, adjusting model transmission parameters with new information would again be a means by which to use our method for proactive intervention planning, followed by updates and adaptation in the event of an outbreak. For example, our post hoc random forest analysis provided some evidence of weak transmission from regressive individuals, in contrast to FeLV dynamics in domestic cats (75). Proactive management planning for FeLV in panthers should, therefore, factor in the risk of transmission from regressively infected individuals, and in the event of an outbreak, update this assumption and management response as new information becomes available about the risk of transmission from regressives.

The suite of tools for inferring transmission networks from infectious agent genomes is rapidly expanding (24). In this study, we used the program Phyloscanner as it maximized the information from our deep sequencing viral data. However, our FIV sequences were generated within a tiled amplicon framework (36, 75), which biases intrahost diversity and limits viral haplotypes (76). Phyloscanner was originally designed to analyze RNA from virions and not proviral DNA, as we have done here. We have attempted to mitigate the effects of these limitations by analyzing several different Phyloscanner outputs to confirm consistency in our results, and by using only binary networks to avoid putting undue emphasis on transmission network edge probabilities, as these are likely uncertain. Further, our primary conclusions from the transmission networks—that age and pairwise distance are important for transmission—are biologically plausible and supported by other literature, as discussed above. Nevertheless, future work should evaluate additional or alternative transmission network inference platforms. In addition, our tiled amplicon framework was not well suited to detection of FIV super- or coinfections, which have been shown to occur in felids (36, 77, 78). Future work with amplification and sequencing methods that are more suited to detection of multiple apathogenic variants could improve confidence for predictions of pathogenic agents and give more insight into the complexities of transmission dynamics.

In addition, ERGMs assume the presence of the “full network” and it is as yet unclear how missing data may affect transmission inferences (39). ERGMs are also prone to degeneracy with increased complexity and do not easily capture uncertainty in transmission events, as most weighted network ERGM (or generalized ERGM) approaches have been tailored for count data [e.g., (79)]. ERGMs may therefore not be the ideal solution for identifying drivers of transmission networks in all systems. Alternatives may include advancing dyad-based modeling strategies (80), which may more easily manage weighted networks and instances of missing data.

Conclusions

Here, we integrated genomic and network approaches to identify drivers of FIV transmission in the Florida panther. This apathogenic agent acted as a marker of close, direct contact transmission, and drivers of FIV transmission were subsequently relevant for predicting the observed transmission dynamics of the related pathogen, FeLV. Further testing of apathogenic agents as markers of transmission and their ability to predict transmission of related pathogens is needed, but they hold promise as a novel tool for proactive epidemic management across host-pathogen systems.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: Data and R code for replication of simulations and analyses is archived and available on Zenodo (https://doi.org/10.5281/zenodo.7025742).

Ethics statement

Ethical review and approval was not required for the animal study because this project did not involve new data collection from vertebrate animals.

Author contributions

MG performed research, analyzed data, and led writing the paper. NF-J, JM, RG, JL, SKr, SKe, RP, and EC performed research. All authors contributed to the design of the research and to writing the paper.

Funding

This research was supported by the National Science Foundation (DEB-1413925, 1654609, and 2030509). MG was supported by the Office of the Director, National Institutes of Health (NIH T32OD010993), the University of Minnesota Informatics Institute MnDRIVE program, and the Van Sloun Foundation. JM was supported by the ACVP/STP Coalition for Veterinary Pathology Fellows and the Linda Munson Fellowship for Wildlife Pathology Research. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.