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.

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).

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.

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)].

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.

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)].

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].

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).

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).

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_k where k is number of nearest neighbors considered: T₃ = 20, p = 0.049; T₅ = 32, p = 0.028; T₇ = 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).

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).

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.

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.

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.