Virus‐Dependent Geographic Structure of Co‐Circulating Viruses in a Single Bat Species

TL;DR

This study shows that different viruses in the same bat species spread differently based on their infection strategies and human activity, not just bat movement.

Contribution

The study reveals that virus-specific infection strategies and human factors, rather than host ecology alone, shape viral population structures.

Findings

Betaherpesvirus and rabies virus spread is influenced by bat travel distance.

Dependoparvovirus is constrained by human-related factors.

Human travel difficulty affects all viruses but with varying impacts.

Abstract

Understanding the spatial spread of viruses within wildlife populations is often a key component of disease management efforts. Viral spread is likely constrained by host ecology, but inter‐virus differences in infection strategy might allow some viruses to overcome these constraints, leading to divergent population structures within a common host environment. We studied the phylogeographic structure of six virus taxa (dependoparvovirus, deltavirus, mastadenovirus, betaherpesvirus and two lineages of rabies virus) circulating in common vampire bats (
Desmodus rotundus
) in Peru, finding that viral population structure is inconsistently constrained by host ecology. Specifically, while bat travel distance structured the genetic diversity of betaherpesvirus and two lineages of rabies virus, other viruses were instead constrained by anthropogenic factors (dependoparvovirus), or had weakly defined population structure (mastadenovirus). The genetic structure of all viruses was affected by a measure of human travel difficulty between sites, but effects varied in size and direction. Distinct drivers of viral population structure within the same host species imply that virus infection strategy can outweigh host ecological connectivity, acting as a key determinant of geographic spread. Because barriers to gene flow generalise poorly between viruses, whether a tractable virus can illuminate host population structure or predict the spread of high‐impact viruses depends on individual virus biology.

Figures

Geographic distribution and prevalence of six vampire bat‐associated viruses in Peru. Black points indicate bat roost sites; bar charts indicate the number of positive samples (first number) out of the total screened for each virus from each roost site (second number), with the site name below each chart. Numbers of sequences from some viruses, particularly betaherpesvirus, are higher than the number of positive samples due to multiple strain infections. Blue and purple points indicate rabies lineage 1 and 4 livestock detection sites, respectively. Black outlines show administrative departmental borders. DrAdV =  D. rotundus mastadenovirus. DrBHV =  D. rotundus betaherpesvirus. DrDV =  D. rotundus deltavirus. DrPV =  D. rotundus dependoparvovirus. DrRV1 =  Desmodus rotundus rabies virus lineage 1. DrRV4 =  D. rotundus rabies virus lineage 4.

Phylogeographic structure of six vampire bat associated viruses. Tip colours indicate ecological region. Scale bars indicate number of nucleotide substitutions. Nodes with at least 75% bootstrap support are marked with a black dot for legibility; Trees with complete bootstrap values are presented in Figures S3–S8. (A) Mastadenovirus tree of 36 partial hexon protein sequences of 478 nt. (B) Betaherpesvirus tree of 323 UL55 glycoprotein B sequences of 150 nt. (C) Deltavirus tree of 43 partial delta antigen sequences of 192 nt. (D) Dependoparvovirus tree of 83 sequences of 951 nt. (E) Rabies virus lineage 1 tree of 86 partial nucleoprotein sequences of 1350 nt. (F) Rabies virus lineage 4 tree of 171 partial nucleoprotein sequences of 1350 nt.

Different variables explain the phylogeographic structure of distinct vampire bat‐associated viruses. Each panel shows the effect sizes of each variable for each virus, with the mean estimate as a point, the 80% highest posterior density estimate (HPDI) as a thick line, and the 95% HPDI as a thin line. Estimates are coloured grey when the 80% HPDI overlaps with 0, medium teal when the 80% HPDI does not overlap with 0 but the 95% HPDI does, and dark teal when the 95% HPDI does not overlap with 0. Neither rabies lineage could be associated with mitochondrial data.

Tables

TABLE 1. Viruses included in analysis show a diverse range of genomes and life history strategies.

Virus	Genome	Samples	Biologically relevant features a
D. rotundus rabies virus, lineage 1	−ssRNA	86	Highly pathogenic, saliva transmitted, found in many departments of Peru (Streicker et al. 2016).
D. rotundus rabies virus, lineage 4	−ssRNA	171	Highly pathogenic, saliva transmitted, found mostly in the Southern Andes of Peru (Cárdenas‐Canales et al. 2022).
D. rotundus deltavirus	−ssRNA	42	Presumed satellite virus requiring a yet unknown helper virus to complete infection cycle (Bergner, Orton, et al. 2021).
D. rotundus Betaherpesvirus	dsDNA	323	Minimally pathogenic with lifelong infections characterised by latency and reactivation (Griffiths et al. 2022).
D. rotundus mastadenovirus	dsDNA	36	Likely host specific, but some cases of cross‐species transmission recorded (Hackenbrack et al. 2017; Kohl et al. 2012). Possible enteric infection in bats (Kohl et al. 2012).
D. rotundus dependoparvovirus	ssDNA	83	Minimally‐pathogenic satellite virus often requiring mastadenovirus or herpesvirus helper to replicate (Flotte and Berns 2005; Lau et al. 2017). Often found in alimentary tract samples (Zhu et al. 2020), so might be an enteric infection.

TABLE 2. Mantel tests correlating the mean branch lengths between viruses at each site and each explanatory variable.

Parameter		DrAdV	DrBHV	DrDV	DrPV	DrRV1	DrRV4
Human accessibility	p	0.353	0.039	0.035	0.069	0.306	0.618
	r	0.149	0.214	0.313	0.181	0.063	0.015
Bat travel distance	p	0.574	< 0.001	0.595	0.025	0.171	0.749
	r	−0.035	0.748	−0.091	0.207	0.097	0.009
PC1	p	0.141	0.537	0.284	0.220	0.377	0.805
	r	0.324	−0.018	0.196	0.059	0.015	0.003
PC2	p	0.574	0.047	0.284	0.017	0.171	0.964
	r	−0.026	0.210	0.233	0.343	0.011	−0.044
Mitochondrial distance	p	0.366	0.169	< 0.001	0.049	NA	NA
	r	0.137	0.108	0.695	0.226	NA	NA

Full text

Introduction

1

Understanding the dynamics and drivers of virus transmission is essential for predicting disease risk and guiding effective action to minimise spread and spillover (Sokolow et al. 2019; Olival et al. 2020; Meza et al. 2022). These dynamics can be affected by both host (e.g., behaviour, dispersal, reproductive strategy) and viral ecology (e.g., host range, transmission by arthropod vectors, acute vs. persistent infection) (Jolles et al. 2021; Wille et al. 2024). Differences in the transmission dynamics of rapidly evolving viruses can manifest as distinct patterns in phylogenetic analysis (Grenfell et al. 2004). For example, acute respiratory viruses such as human influenza virus and SARS‐CoV‐2 show source‐sink and wave‐like spread driven by travel and human density, whereas sexually transmitted or blood‐borne chronic viruses such as HIV and hepatitis C tend to be genetically compartmentalised within certain geographic areas and risk groups (Takebe et al. 2004; Pybus et al. 2005; Holmes 2009; Li et al. 2021). Factors other than transmission mode, such as variation in the average age of infection, viral host range (including transmission by arthropod vectors), and environmental transmission may introduce further unpredictability in viral phylogeographic structure (Salemi et al. 1999; Sejvar 2003; Leggett et al. 2013; Leung 2021). The presence of co‐circulating viruses may also have direct or indirect effects on viral spread through cross‐immunity or direct interference (Rohani et al. 2003; Nickbakhsh et al. 2019). For example, ‘satellite’ viruses (viruses that require a co‐infecting helper virus to complete their life cycle) are established and emerging human health threats and may represent an extreme case (Netter et al. 2021; Ho et al. 2023). Compared to viruses that do not obligately require co‐infection, satellite virus spread may be more influenced by the distribution, prevalence and host range of the helper virus than the host.

To date, phylogenetic comparisons of viruses within common host environments have largely focused on human pathogens. While insightful for understanding human risk, distinct features of human ecology mean the spread of human viruses are unlikely to reveal drivers of viral spread in natural populations. For example, human virus spread is influenced by hyper‐connected populations supported by air travel, which will not generalise to viruses of most other species. Targeted collection of viruses from clinical settings and medical (e.g., vaccination) and non‐medical (e.g., quarantine) interventions further introduce sampling biases and artificial constraints on transmission, which are rare in natural populations. Relative to humans, the higher sensitivity of wildlife to landscape barriers, the rarity of interventions to disrupt natural viral spread, and the more integrated multi‐species systems make wildlife viruses well placed to illuminate the relative roles of environment and viral life history on spatial infection dynamics in nature.

Wildlife viruses are also an increasing concern for public and veterinary health, making their spatial spread important in its own right (Jones et al. 2008; Sims and Brown 2016; Jolles et al. 2021; Li et al. 2023; De Campos et al. 2024). Yet, because many viruses are difficult to sample, more easily studied wildlife viruses have been proposed to serve as markers of host population demographics and connectivity (Biek et al. 2006; Fountain‐Jones et al. 2017). Similarly, host population structure and connectivity have been proposed as informative factors to anticipate the spread of both endemic and invading viruses (Olival et al. 2020; Streicker et al. 2016). However, the extent that phylogeographic patterns generalise across virus species in the same host remains a key unknown which will govern the appropriateness of using viruses as proxies for host or other virus dispersal (Carver and Lunn 2020). Specifically, if biological differences among viruses erase the imprint of host species behaviour, then viral population structures may diverge significantly from the host and from other viruses, limiting their value as proxies. Detecting divergent population structure among viruses requires investigating multiple viruses within a shared host species, environment and temporal window.

Common vampire bats ( Desmodus rotundus ) constitute a tractable system to investigate the population structures of co‐circulating viruses. Vampire bats occur across most of Central America, South America, and southern North America in diverse habitats including the rainforests, mountains, and coastal deserts (Greenhall et al. 1983). Obligate blood feeding puts vampire bats into regular and direct contact with humans, domestic animals and wildlife, facilitating cross‐species transmission. Indeed, in many parts of Latin America, vampire bats are the primary source of human rabies (Schneider et al. 2009). Vampire bats also host numerous non‐rabies viruses, some of zoonotic concern (Bergner et al. 2020; Bergner, Orton, et al. 2021; Griffiths et al. 2022). Phylogeographic studies to date have shown isolating effects of mountains on both vampire bats and some of their viruses (Streicker et al. 2016; Griffiths et al. 2020). To assess the cross‐viral predictive value of a given host and virus genetic dataset, we need to understand whether barriers to host gene flow constrain all co‐circulating viruses, or whether some viruses escape these barriers. Currently, this is unknown.

Here, we compare the phylogeographic structure of six viruses which co‐circulate in the same populations of common vampire bats in Peru, which contains significant geographical, ecological and prey species diversity, as well as varying levels of urban and agricultural impact on the landscape. The viruses we investigate have a range of life histories that may cause differences in their transmission dynamics and therefore population structure, including: (1, 2) two strains of a directly transmitted RNA virus which causes acute lethal infections (rabies virus) (Cárdenas‐Canales et al. 2022); (3) a presumed RNA satellite virus with unknown infection strategy, infection duration, and host range ( D. rotundus deltavirus) (Bergner, Orton, et al. 2021); (4) a saliva‐transmitted host‐specific DNA virus with minimal pathogenicity and lifelong within‐host persistence ( D. rotundus betaherpesvirus) (Griffiths et al. 2023); (5) a double‐stranded DNA virus with unknown ecology but potential direct transmission and varying host range ( D. rotundus mastadenovirus) (Kohl et al. 2012); and (6) a single‐stranded DNA presumed satellite virus with unknown ecology and host range but potential direct transmission from enteric infections ( D. rotundus dependoparvovirus) (Zhu et al. 2020) (Table 1). We carried out virus‐specific modelling of molecular sequence data to determine whether these differences in infection strategies affect viral population structure. We hypothesised that, if host connectivity primarily governs viral population structure, a common set of predictors of viral genetic structure (e.g., geographic or host genetic distance) would emerge, whereas stronger effects of virus life history would cause contrasting patterns of viral population structure and different underlying predictors.

Methods

2

We selected six viral taxa that, based on previous metagenomic sequencing or surveillance studies, were suspected to be relatively common in bats (or for rabies, livestock, since rabies virus sequences are rare in sampled bats) and were widely distributed in Peru. Sequences from three viruses (deltavirus, rabies virus and betaherpesvirus) were assembled from previous studies (Bergner, Orton, et al. 2021; Griffiths et al. 2022; Streicker et al. 2016). Sequences from mastadenovirus and dependoparvovirus were generated by nested polymerase chain reaction (PCR) and sequencing (see below). The full dataset included 741 sequences from six viruses across 30 bat roosts and over 200 unique rabies cases. Altogether, data were collected from 13 departments and all major ecological regions of Peru (coastal desert, mountains, and tropical rainforest), although not every virus was represented in every department or ecological region (Figure 1; Figure S1). Rabies sample collection dates ranged from 1997 to 2012, and non‐rabies samples were collected between 2009 and 2018, with most collected between 2015 and 2018 (see below for details).

Geographic distribution and prevalence of six vampire bat‐associated viruses in Peru. Black points indicate bat roost sites; bar charts indicate the number of positive samples (first number) out of the total screened for each virus from each roost site (second number), with the site name below each chart. Numbers of sequences from some viruses, particularly betaherpesvirus, are higher than the number of positive samples due to multiple strain infections. Blue and purple points indicate rabies lineage 1 and 4 livestock detection sites, respectively. Black outlines show administrative departmental borders. DrAdV =  D. rotundus mastadenovirus. DrBHV =  D. rotundus betaherpesvirus. DrDV =  D. rotundus deltavirus. DrPV =  D. rotundus dependoparvovirus. DrRV1 =  Desmodus rotundus rabies virus lineage 1. DrRV4 =  D. rotundus rabies virus lineage 4.

Previously Published Viral Sequence Datasets

2.1

The D. rotundus deltavirus dataset consisted of 43 partial delta antigen sequences (192 nucleotide, nt), and came from saliva and blood samples collected between 2016 and 2017 from 10 sites in three departments of Peru: Ayacucho, Lima and Cajamarca (Bergner, Orton, et al. 2021). The D. rotundus betaherpesvirus dataset included 323 partial UL55 glycoprotein B sequences that were generated by amplicon sequencing saliva samples from 23 vampire bat roosts (sampled between 2009 and 2018) on an Illumina sequencing platform. The high frequency of co‐infection in this dataset allowed separation of multiple viral strains within the same sample but restricted sequences to 150 nt reads (Griffiths et al. 2022). Both the deltavirus and betaherpesvirus data were limited to short sequence lengths due to practical constraints: the deltavirus fragment was short due to challenging PCR design and concerns over false negatives with a longer sequence (Bergner, Orton, et al. 2021), and the betaherpesvirus fragment was short due to the challenges of building long non‐chimeric consensus sequences from short read data in individuals infected with multiple strains (Griffiths et al. 2022). The rabies virus dataset included 264 partial nucleoprotein sequences (1350 nt) collected from livestock brains in 13 departments of Peru (Streicker et al. 2016). These sequences represent two main clades of vampire bat‐transmitted rabies virus, previously described as lineage 1 (86 sequences collected between 1997 and 2012) and lineage 4 (171 sequences collected between 1997 and 2012). Although not collected directly from bats, livestock do not transmit rabies virus, so each sequence represents a single transmission event from bats, and these livestock sequences are therefore considered proxies for rabies virus circulating in the local D. rotundus populations (Streicker et al. 2016). Methods on sampling, sample processing and sequencing are provided in the original articles describing the samples analysed here (Bergner, Orton, et al. 2021; Bergner et al. 2020; Griffiths et al. 2022; Streicker et al. 2016).

Amplification and Sequencing of Additional Vampire Bat‐Associated Viruses

2.2

We selected viruses which were found in at least six out of 31 pools of existing metagenomic sequencing data and spread across at least three administrative departments of Peru (Bergner, Mollentze, et al. 2021). Viruses known to have multi‐strain infections or extremely high diversity were generally excluded since Sanger sequencing would have created chimeric sequences due to multiple strains being present in the PCR product, making phylogeographic analysis challenging or impossible (Moreira Marrero et al. 2021). This led to the exclusion of gammaherpesvirus and picornavirus. Exceptions were allowed for Desmodus rotundus betaherpesvirus, since the dataset for this virus was generated by next‐generation sequencing, allowing us to separately analyse co‐infecting strains (see above), and for a small number of mastadenovirus and dependoparvovirus infections, where we were able to use cloning to separate co‐infections (see below). In total, our additional PCR screening allowed us to add a vampire bat‐associated mastadenovirus and a vampire bat‐associated dependoparvovirus to the previously published datasets (Table 1).

Mastadenovirus and dependoparvovirus PCR primers were designed by aligning contigs from metagenomic sequencing of vampire bat rectal swabs (Bergner, Mollentze, et al. 2021; Bergner et al. 2020) to bat mastadenovirus A (accession: LC385828.1) and bat dependoparvovirus (accession: GU226971) genomes using Geneious Prime 2021.2.2 (http://www.geneious.com). Consensus sequences were created for the mastadenovirus hexon and dependoparvovirus VP2 proteins, and a range of primers were designed using Geneious and Eurofins Primer Design Tool (eurofinsgenomics.eu). Proteins were chosen by their use in other phylogenetic studies (Clegg et al. 2011; Dhingra et al. 2019) and the presence of high‐quality contigs in metagenomic data for primer design. Primers and cycling conditions are available in Table S1.

For both new viruses, a total of 289 D. rotundus nucleic acid extracts from rectal swab samples collected between 2009 and 2017 were screened using nested PCR, including all previously published samples and 52 unpublished samples collected as part of the same project (Bergner et al. 2020) (Table S1). PCR products were analysed using agarose gel electrophoresis, and positive products were Sanger sequenced (Eurofins Scientific). Sequences were aligned and manually trimmed using Geneious. Candidate viral sequences were confirmed using nucleotide BLAST analysis against the NT database (Altschul et al. 1990). One mastadenovirus PCR product and three dependoparvovirus PCR products were excluded from phylogenetic analysis due to poor sequence quality. One mastadenovirus and five dependoparvovirus sequences appeared to be mixed infections based on double peaks in chromatograms at the same nt sites in both the forward and reverse sequencing directions. These apparently co‐infected samples were cloned into plasmids using the TOPO cloning kit (Invitrogen) following the manufacturer's recommended conditions, and five cloned plasmids from randomly chosen individual bacterial colonies of each original PCR product were sequenced and processed as described above. A preliminary phylogenetic analysis showed that the sequences from the cloned plasmids often formed two to three closely related clusters in different parts of the phylogenetic trees, each of which we assumed reflected a single infection from that strain. Subsequent analyses therefore included a single randomly chosen sequence from each cluster to provide only one sequence per infecting strain.

Amplification and Sequencing of

Desmodus rotundus Cytochrome B Gene

2.3

We used vampire bat population mitochondrial structure as a proxy for host connectivity. Complete mitochondrial cytochrome B (cytB) sequences were generated for several sites for which data was not already available, targeting 15 individuals per site (aside from LR2, where only eight samples were available). Following previously described protocols (Martins et al. 2007; Streicker et al. 2016), DNA from individual bat wing biopsies was extracted using the DNEasy Blood & Tissue Kit (Qiagen), amplified by PCR and Sanger sequenced (Eurofins Scientific). A maximum likelihood phylogenetic tree of all cytB sequences (ranging from 741 to 1140 nt) was made with IQ‐Tree using default settings for automatic evolutionary model selection by Bayesian information criterion (selected model: TN + F + R6) (Minh et al. 2020), and pairwise branch lengths were extracted and averaged between each pair of sites, including within‐site, using ape v.5.7–1 (Paradis and Schliep 2019) in R v.4.3.0 (R Core Team 2023).

Collection of Genetic and Environmental Data for Analysis of Viral Population Structure

2.4

We used Bayesian hierarchical models to relate the phylogenetic branch lengths between individual virus sequences with measures of geographic, ecological, and host genetic distance. For each virus, maximum likelihood phylogenetic trees were constructed with IQ‐Tree with automatic evolutionary model selection by Bayesian information criterion as described above (DrAdV: TNe + I. DrBHV: K2P + G4. DrDV: TN + F + I. DrPV: K3Pu + F + G4. DrRV1: TN + F + R3. DrRV4: TN + F + R2) and ultra‐fast bootstrap for 1000 repetitions (Minh et al. 2020). Pairwise branch lengths were extracted as above, and trees were visualised using ggtree v.3.9.0 (Yu 2020).

We included three major categories of explanatory variables to explain viral population structure. First, we calculated the connectivity between sampling sites for both bats and humans, as the latter may influence virus transmission as hosts or by translocating infected animals. The host connectivity between each pair of sites was calculated using both maternal population structure (using the cytB data described above) and by a least‐cost paths model previously developed for this environment (Benavides et al. 2016), where travel costs increase linearly with elevation (hereafter ‘bat travel distance’). Elevations > 4000 m were deemed impassable due to physiological constraints based on to the distribution of bat roosts observed and the success of this threshold in other models in this system (Benavides et al. 2016; Greenhall et al. 1983; Ribeiro et al. 2023). Although elevation is not the only variable that might affect vampire bat movement, other variables, such as mean forest cover and ruggedness of terrain, were included separately to account for environmental differences (described below), allowing us to disentangle any effects of climate from distance between sites and providing an easily interpretable measurement of travel cost. As bat travel distance is correlated with vampire bat relatedness using bi‐parentally inherited nuclear markers (r = 0.69) and nuclear and mitochondrial markers weakly correspond in this species, we considered bat travel distance as a proxy for male and female dispersal (Martins et al. 2007; Streicker et al. 2016; Griffiths et al. 2022). This distance calculation was applied identically to both bat roost sites (non‐rabies viruses) and livestock sampling sites (rabies viruses). We consider livestock samples to be proxies for bat roost sites due to the short foraging range of vampire bats (< 5 km) (Greenhall et al. 1983; Ribeiro et al. 2025) and the infrequency of long‐distance livestock movement in Peru. Human travel accessibility between each site was estimated by the same least‐cost paths method using the Malaria Atlas Project friction surface (Weiss et al. 2020), which quantifies human travel effort using features such as roads, water crossings and elevation (hereafter ‘human accessibility’).

Second, we gathered factors describing the environments where viruses were found, reasoning that precipitation, temperature and plant cover and type could affect gene flow if viral persistence in the environment is important for transmission (such as higher humidity facilitating transmission [Lowen et al. 2007]), if certain alternate hosts are present only under certain ecological conditions, or if environmentally driven behavioural differences affect viral transmission (Herrera and Nunn 2019; Lowen et al. 2007) (Table S2). In addition, human population density and footprint (i.e., scale of impact humans have had on the land) and livestock density were recorded, as these may be important if a virus uses humans or domestic animals as alternate hosts, or if these organisms cause ecological changes that alter bat dispersal and thus virus transmission (Viana et al. 2023) (Table S2). Data from all available livestock species were included in this analysis (Table S2), as D. rotundus is a generalist feeder (Brown and Escobar 2023).

Third, we calculated the number of years between sampling dates for each pair of sequences to account for measurable viral evolution during the study period, such as that observed in the rabies lineages. Although the sequences for other viruses were collected on a much shorter time scale, making measurable evolution unlikely, the difference in years was included as an explanatory variable in all virus models for consistency and to account for potential differences in the spatial allocation of sampling effort across years (Table S2). Additionally, although both lineages of rabies virus were sampled over a much larger time span (1997–2021 vs. 2009–2018), the time of sampling is unlikely to otherwise bias the analyses, as bat demographics are not expected to have changed substantially over these time periods. The environmental variables are also estimates from 1 year each, making them similarly resolved for all viruses.

When appropriate for our hypotheses, environmental variables of different original spatial resolutions were rescaled using bilinear interpolation (see Table S2 for original resolutions). Specifically, variables associated with bat foraging (livestock and human density) were resampled to 5 km^2^ based on knowledge of typical vampire bat foraging distances (Greenhall et al. 1983). Variables associated with the natural environment (e.g., slope, all climatic variables) were resampled to 0.1 km^2^ (the finest resolution present) to reflect the hypothesised influence of local conditions. The mean value of each environmental variable in a 10 km buffer around the site (based on the maximum bat foraging distance [Greenhall et al. 1983]) was then extracted (Table S2). Although resampling can increase spatial autocorrelation, we reduced this issue by averaging values from the 10 km buffer around each sampling site and using widely spaced sites (over 5 km apart with few exceptions). Environmental variables were normalised and subjected to principal components analysis (PCA). Subsequent analyses relating viral genetic divergence to environmental characteristics used the absolute value of the difference in values of the first two principal components (PCs; together capturing 70% of the observed variation) between each pair of sites as explanatory variables. PC1 contained bioclimatic variables including altitude, temperature and precipitation, as well as horse density. PC2 contained some livestock density variables (specifically of chickens and pigs) as well as human‐related variables such as population density and human footprint (details of PCA available in Table S4 and Figure S2). An effect of PC1 could therefore indicate preferential viral dispersal between similar biomes due to host preferences or local bioclimatic conditions that facilitate viral circulation. Effects of PC2 could indicate that humans or chicken and pigs, or the human movement of these animals, facilitate transmission between similar, human‐dominated environments which are otherwise isolated by bat movement alone. The potential interpretations of all variables included in this analysis are summarised in Table S3.

Statistical Modelling

2.5

Infection prevalence for newly studied viruses was calculated with 95% binomial confidence interval using R v.4.3.0 (R Core Team 2023). Binomial generalised linear models (GLMs) were fitted to explain the detection of mastadenovirus and dependoparvovirus in each sample, with bat age, sex and reproductive status as explanatory variables, and site as a random effect. Bat age was coded as juvenile or adult based on whether the epiphyseal and diaphyseal bones in the wing were fused (Anthony 1988). Reproductive status was coded as binary, with females coded as reproductive if they were pregnant or lactating, and males coded as reproductive if they had enlarged testicles (Streicker et al. 2012). Because mastadenoviruses are hypothesised candidate helper viruses for dependoparvoviruses, we also included the presence or absence of the other virus as an additional binary variable. Model selection was performed using the corrected Akaike's Information Criterion (AICc) (Akaike 1973), with ∆AICc > 2 indicating differences among candidate models. Fixed effect selection was performed using MuMIn v.1.47.5 (Barton 2023), with all possible combinations of explanatory variables tested and candidate models ranked by AICc, with the best performing model presented.

Preliminary analyses of the relationships between each explanatory variable and mean viral divergence between pairs of sites used Mantel tests with 9999 permutations within the vegan v.2.6–4 package (Oksanen et al. 2022). Mantel test p‐values were adjusted for multiple testing using the Benjamini–Hochberg correction (Benjamini and Hochberg 1995). Because the Mantel tests do not account for phylogenetic non‐independence, we used normally distributed Bayesian hierarchical models to explore relationships between the environmental covariates described above and the pairwise phylogenetic distances among viruses from different sites and identify potential drivers of spatial spread. All explanatory variables were normalised before being included in the analysis. Random effects were site pair (the combination of each pair of sites where viruses were sampled, accounting for the fact that the branch lengths between these sites will have identical explanatory variables), SiteA and SiteB (the sites each virus was sampled in, accounting for any particularly divergent sites), SeqA and SeqB (the sequences to and from which the branch length was calculated, accounting for the fact that some individual sequences may be particularly divergent), and SeqA and SeqB with a covariance structure calculated from the phylogenetic tree to account for the phylogenetic non‐independence of the sequences. Including both the simple and the phylogenetically corrected random effect structure is intended to account for non‐independence of sequences within a tree (Cinar et al. 2022; Schmidt et al. 2021). Therefore, most model structures were: BranchLength ~ BatTravel + CytB + Accessibility + PC1 + PC2 + (1|SitePair) + (1|gr(SeqA, cov = Tree)) + (1|gr(SeqB, cov = Tree)) + (1|SeqA_2) + (1|SeqB_2) + (1|SiteA) + (1|SiteB).

Slight variations in modelling strategy were required due to idiosyncrasies across the datasets for each virus. Mastadenovirus and dependoparvovirus phylogenies both showed one large and one small clade, and in both cases, analyses were conducted using only sequences from the large clade. Rabies lineages 1 and 4 were analysed independently, and one sequence from each lineage from sites that were implausible to be inhabited by vampire bats (over 4000 m elevation) was excluded, as these are likely the result of seasonal translocations of cattle that were presumably infected at lower elevations. Finally, for rabies virus models, CytB, SiteA, SiteB and SitePair were excluded since rabies sequences originated from livestock and could not be associated with any specific bat roost.

All phylogenetic models were fit with Stan v.2.32.2 in the brms package (Bürkner 2021; R Core Team 2023; Stan Development Team 2023). Models used four chains of 10,000 iterations, with the first 5000 of those iterations discarded as a warm‐up. Model fit was assessed by confirming that Rhat values were 1.01 or below, bulk and tail ESS values were above 400, and there were no divergent transitions (Stan Development Team 2023).

Collinearity in posterior estimates was assessed by examining the pairs plot for each model (Figures S9–S14), but variables where correlation was visible were not removed because they did not lead to poor model performance. Results were presented with the posterior mean estimate alongside the 95% and 80% highest posterior density intervals (HPDI). We considered effect sizes with 95% HPDIs that did not overlap with zero to indicate significant effects with strong confidence, while non‐overlap of the 80% HPDI indicated non‐zero effects with low confidence.

Results

3

Mastadenoviruses and Dependoparvoviruses Are Prevalent and Mutually Associated

3.1

Of 289 D. rotundus rectal swab samples, 37 were positive for mastadenovirus (12.8%, 95% CI: 0.094, 0.171) and 80 were positive for dependoparvovirus (27.7%, 95% CI: 0.228, 0.331). Binomial GLMs showed that the most important predictor for each virus was the presence of the other (mastadenovirus: odds ratio (OR) = 5.78, 95% CI = 2.47–14.37, R ^2^ = 0.163; dependoparvovirus: OR = 1.8, 95% CI = 0.97–2.76, full model R ^2^ = 0.33). The model for dependoparvovirus also included an effect of age, though confidence intervals overlapped with 1 (Tables S5–S7).

Idiosyncratic Population Structure of Co‐Circulating Vampire Bat Viruses

3.2

Phylogenetic trees of most vampire bat viruses indicated spatial structure between the coast, eastern Andean slopes, lowland Amazon, and southern Andean valleys of Peru (Figure 2). For example, both lineages of rabies virus were restricted to areas east of the Andes. Rabies virus lineage 1 was found mostly in the eastern Andean slopes, with only a few sequences in the lowland Amazon and southern Andes, and rabies lineage 4 was restricted to the southern Andes (Figure 2E,F; Figures S3–S8). The deltavirus and dependoparvovirus trees were similarly structured, with many clades containing sequences from the same regions (Figure 2C,D). However, some dependoparvovirus clades contained sequences from different regions, suggesting at least four trans‐Andean dispersal events (Figure 2D). The betaherpesvirus tree contained some clades apparently isolated in the southern Andes as well as clades distributed across the eastern Andean slopes and coastal regions again supporting trans‐Andean viral gene flow (Figure 2B). The mastadenovirus tree showed relatively little spatial structure, with almost all clades containing sequences from multiple regions (Figure 2A). The trees also varied substantially in depth, with some viruses, such as betaherpesvirus and dependoparvovirus, showing more genetic change within the tree than others.

Phylogeographic structure of six vampire bat associated viruses. Tip colours indicate ecological region. Scale bars indicate number of nucleotide substitutions. Nodes with at least 75% bootstrap support are marked with a black dot for legibility; Trees with complete bootstrap values are presented in Figures S3–S8. (A) Mastadenovirus tree of 36 partial hexon protein sequences of 478 nt. (B) Betaherpesvirus tree of 323 UL55 glycoprotein B sequences of 150 nt. (C) Deltavirus tree of 43 partial delta antigen sequences of 192 nt. (D) Dependoparvovirus tree of 83 sequences of 951 nt. (E) Rabies virus lineage 1 tree of 86 partial nucleoprotein sequences of 1350 nt. (F) Rabies virus lineage 4 tree of 171 partial nucleoprotein sequences of 1350 nt.

Correlations of Virus Branch Lengths With Ecological Variables

3.3

Mantel tests correlating average viral branch lengths with each explanatory variable suggested differences in the drivers of viral spatial structure across viruses (Table 2). For example, deltavirus branch lengths were associated with mitochondrial distance but not bat travel distance, and betaherpesvirus branch lengths were strongly associated with bat travel distance (Table 2). Neither lineage of rabies virus showed any relationship between branch lengths and explanatory variables in the Mantel tests.

Ecological Drivers of Population Structure Differ Across Viruses

3.4

The phylogenetic distance models for each virus confirmed that different predictors influenced spatial genetic structure across viruses (Figure 3, Table S8). Human accessibility had a non‐zero effect on every virus, although this effect was less confident in mastadenovirus, with the 95% HPDI still overlapping with zero (Figure 3, Table S8). However, the effects of accessibility were not all in the same direction: effects were positive for betaherpesvirus, deltavirus, dependoparvovirus, and rabies virus lineage 4, but negative for mastadenovirus and rabies virus lineage 1 (Table S8, Figure 3). In contrast, bat travel distance had a strong and high‐confidence positive effect on both rabies virus lineages and betaherpesvirus, a relatively small negative effect on deltavirus, and had no effect on other viruses (Table S8, Figure 3). Another common predictor was the difference in PC2, which had a positive effect on both dependoparvovirus and rabies virus lineage 4, but a small negative effect on betaherpesvirus and rabies virus lineage 1 (Table S8, Figure 3). In comparison, the effect of the difference in PC1 was only non‐zero for mastadenovirus and dependoparvovirus, and then only with low confidence (Table S8, Figure 3). Difference in year had a positive effect on mastadenovirus, both rabies lineages (which are known to be measurably evolving [Streicker et al. 2016]), and a relatively weak and low‐confidence effect on dependoparvovirus (Table S8, Figure 3). Mitochondrial distance did not confidently influence viral genetic structure in any model, but it had a weak negative effect on dependoparvovirus and a large but highly uncertain effect on deltavirus (Table S8, Figure 3). The pairs plots for all models indicated that no parameter estimates were correlated with each other except for bat travel distance and human accessibility in both rabies lineages (Figures S9–S14).

Different variables explain the phylogeographic structure of distinct vampire bat‐associated viruses. Each panel shows the effect sizes of each variable for each virus, with the mean estimate as a point, the 80% highest posterior density estimate (HPDI) as a thick line, and the 95% HPDI as a thin line. Estimates are coloured grey when the 80% HPDI overlaps with 0, medium teal when the 80% HPDI does not overlap with 0 but the 95% HPDI does, and dark teal when the 95% HPDI does not overlap with 0. Neither rabies lineage could be associated with mitochondrial data.

Discussion

4

Given that viruses rely on their hosts for transmission, host dispersal and population connectivity necessarily influence viral spatial spread. Many studies have exploited this expected correspondence to make biological inferences that would be impractical using other methods. For example, fast evolving, host‐associated feline retroviruses have been used to infer aspects of wild felid demographic history, behaviour, and spatial connectivity which could not have been easily inferred from either host genetics or field studies (Biek et al. 2006; Fountain‐Jones et al. 2017). In contrast, host genetic connectivity has been used to infer the spatial spread of bat‐associated rabies viruses to currently rabies‐free areas (Streicker et al. 2016). However, variation in transmission route and host specificity across viruses means that the extent to which host dispersal governs the spread of different viruses remains unclear (Carver and Lunn 2020). Analysing multiple viruses in the same host species provides a rare opportunity to explore whether predictors of viral gene flow are consistent across viruses. Specifically, we demonstrated idiosyncratic phylogeographic structures and ecological drivers of gene flow in six viruses across the same vampire bat host populations. Our results suggest that differences in viral biology lead to varying dependence on host movement for viral spread at the landscape scale, with practical implications for using viruses as proxies for host movement.

To our knowledge, all of the viruses we studied are maintained in nature by vampire bats, creating an expectation that bat connectivity should predict viral gene flow and spatial structure. Interestingly, proxies for host connectivity shaped the spatial genetic structure of several, but not all, vampire bat‐associated viruses. Specifically, isolation by bat travel distance (a proxy for overall population connectivity) was positively correlated with the genetic structure of betaherpesvirus and both lineages of rabies virus, but had a weak negative effect on deltavirus and no effect on dependoparvovirus or mastadenovirus. We further explored the effects of bat mitochondrial DNA divergence as an alternative metric of bat population structure, but effects were consistently absent or uncertain across viruses. For rabies virus and betaherpesvirus, the stronger influence of bat travel distance over mitochondrial structure may reflect the importance of male bats for the spatial spread of both viruses, as previously argued for rabies virus, and is consistent with previous studies indicating congruence of both viruses with host population structure (Streicker et al. 2016; Griffiths et al. 2020).

One possible exception to the absent effects of mitochondrial population structure was a positive effect of mitochondrial distance on deltavirus in both the Mantel tests and the phylogenetic models, with the latter result highly uncertain. If substantiated, this effect (together with the putatively negative effect of our landscape metric of overall bat travel distance) could point to female‐biased transmission, as might occur if transmission is predominantly vertical (from mothers to offspring). Vertical transmission is presumed rare in human deltavirus due to the typically low viremia of the hepatitis B helper virus (Ramia and Bahakim 1988; Sellier et al. 2018; Aliasi‐Sinai et al. 2023), but it could be enabled in bats through use of an alternative helper virus with a greater predilection for vertical transmission or if bat deltaviruses can spread without helper viruses (Paraskevopoulou et al. 2020; Bergner, Mollentze, et al. 2021).

Although the general pattern which emerged from our analysis was that different factors shaped the spatial genetic structure of different viruses within the same host population, human accessibility had non‐zero effects for every virus (Figure 3). Crucially, however, the direction of these effects varied across viruses. Negative effects of human accessibility, as observed for mastadenovirus and rabies virus lineage 1, suggest a high level of connectivity between viruses in areas that are difficult for humans to travel between, and conversely, high diversity in areas that are easy for humans to travel between. This indicates transmission with mechanisms that are highly divergent from human movement. In contrast, positive effects observed for dependoparvovirus and rabies virus lineage 4 suggest that these viruses are more closely related in areas that are easier for people to travel between, which could suggest human‐mediated transmission or otherwise similar constraints upon both humans and these viruses.

Human accessibility had different effects even between the two lineages of rabies virus, highlighting that even highly related viruses may have different drivers of spatial spread. This difference also suggests that rabies dynamics may vary between the valleys in the Southern Andes (where rabies lineage 4 was sampled), with the virus moving across the landscape in ways very similar to how humans can move, and the rest of the country (where rabies lineage 1 was sampled), with the virus moving somewhat inversely to how humans can move (Streicker et al. 2016). This might occur if rabies‐infected bats disperse in different ways in mountain valleys compared to more lowland or forested areas as a response to the different environments, or if accessibility has different effects at the smaller spatial scales observed in rabies lineage 4 compared to the large distances seen in rabies lineage 1. Although the accessibility and bat travel distance parameter estimates were somewhat correlated in both rabies lineages, potentially due to sampling livestock sites rather than bat roosts, this did not affect model performance, indicating that even when bat travel distance is accounted for, accessibility has a detectable and meaningful effect.

Examining the determinants of the spatial structure of viruses that were independent of proxies for host connectivity may provide new insights into their biology. For dependoparvovirus, viral spatial structure was confidently correlated with human‐associated variables (i.e., human accessibility and PC2), weakly correlated with differences in the environment (PC1), and phylogenetic analysis indicated at least four trans‐Andean viral dispersals, which are unlikely to be a result of natural vampire bat movement. Although human translocations of infected wildlife can facilitate long‐distance viral dispersals, such as the movement of rabies‐infected raccoons on garbage trucks (Wilson et al. 1997), similar translocations are rarely reported for vampire bats and unlikely to cause this effect (Constantine 2003). However, human movements of infected domestic animals might enable viral dispersal. Indeed, dependoparvoviruses have been reported in several mammals and birds, including in species that are regularly fed upon by vampire bats, and appear capable of replicating in the cells of most vertebrates (Flotte and Berns 2005; Pénzes et al. 2015; Lau et al. 2017; Becker et al. 2018; Dai et al. 2022). Although reliance of dependoparvoviruses on either mastadenovirus or herpesvirus as helper viruses adds an additional constraint on inter‐species transmission (Flotte and Berns 2005; Lau et al. 2017), both of these taxa are common in wild and domestic species (Woźniakowski and Samorek‐Salamonowicz 2015; Glavina et al. 2022), and other satellite viruses have broad host ranges despite reliance on helpers (Zarzyńska‐Nowak et al. 2020; Bergner, Orton, et al. 2021). Additionally, scientific understanding of dependoparvoviruses is changing: previously thought to be non‐pathogenic in humans, these viruses are now understood to cause disease in certain circumstances (Ho et al. 2023), suggesting that they may be capable of much more than currently known. Efforts to detect bat‐associated dependoparvoviruses in livestock, such as the pigs and chickens represented by PC2, would be an important next step to understand whether human translocation of livestock enables viral spread between otherwise isolated vampire bat populations.

For mastadenovirus, divergence was only confidently associated with the difference in year, which was surprising given the generally slow evolutionary rates of dsDNA viruses and the short time scale of our sampling (only 2015–2016) (Duffy et al. 2008). Therefore, we suspect this may be an artifact, perhaps attributable to sites sampled in different years containing different strains. Interestingly, a low‐confidence positive effect of PC1 on mastadenovirus was also detected in the phylogenetic distance models and the Mantel tests. If verified with a larger dataset, it could indicate that viruses in similar environments are more similar than viruses close in geographic space but in different ecological spaces. This could occur if these viruses also circulate within other bat species that are not as widespread, but that co‐roost with D. rotundus , which has been observed in some areas (Griffiths et al. 2020), and would be consistent with the low‐confidence negative relationship with human accessibility suggesting that human movement is an unlikely driver for spread in this virus.

The inconsistent effects of host connectivity, accessibility, and environmental conditions across viruses imply that a virus's value as a proxy for host or other virus dynamics varies, potentially due to biological differences. Considering this diversity, we suggest choosing well‐studied viruses with known biology when using them to infer host dynamics, as has been done successfully with big cats and feline immunodeficiency virus (Biek et al. 2006; Fountain‐Jones et al. 2021), and interpreting the results in the context of each virus's unique infection dynamics. Our analysis also identified potential opportunities for extrapolation, particularly between viruses that have shared transmission modes, such as rabies virus and betaherpesvirus. Rabies virus causes an acute, lethal infection, while betaherpesvirus causes a lifelong, benign infection, but both are saliva transmitted (Griffiths et al. 2020), suggesting that this may be a more important factor in determining population structure than other aspects of viral infection strategy, such as pathogenicity. This result is consistent with the literature on human viruses, which show similarities between blood‐borne and respiratory transmission, even when the viruses in question are phylogenetically distant and have diverse pathobiology (Takebe et al. 2004; Pybus et al. 2005; Holmes 2009; Li et al. 2021). In contrast, mastadenovirus and dependoparvovirus had no robust association with bat connectivity or other viruses, suggesting low utility as generalisable markers of bat or virus spread and that factors beyond host movement governed the spread of these viruses. If generalisable patterns can be found between viruses of different transmission modes, then this could lead to these viruses being used to capture and compare different biological processes in their hosts, such as different types of contact.

The comparative phylogeographic approach we adopted has limitations. First, the sample size, sequence lengths and spatial distribution varied between viral taxa since they depended on virus detections from samples collected in a variety of field studies. While we attempted to minimise these effects by setting inclusion criteria for relatively common and widespread viruses, the reliability and confidence in the analyses necessarily varied across viruses. However, despite these differences, we were able to detect effects associated with spatial structure for most viruses, and we have no reason to believe that variation in sequence length or sample size would have led to selection of different ecological drivers of viral spatial structure. Indeed, betaherpesvirus and deltavirus had similar geographic extents of sampling but non‐overlapping predictors of spatial structure. Second, while some viruses had strong bootstrap support in phylogenetic analyses and clear statistical drivers of spatial structure (e.g., dependoparvovirus), others (e.g., betaherpesvirus and deltavirus), relied on short sequences that may have limited the resolution of our phylogenetic analysis and power of our statistical modelling. In contrast, both rabies lineages had a large sample size, sequence length, and more continuous spatial distribution than most other viruses, potentially providing a higher resolution of data that potentially improved our ability to detect finer scale drivers of viral spread (Streicker et al. 2016). However, given the consistency of some of these models with other studies on these viruses (Streicker et al. 2016; Griffiths et al. 2020), we expect that with increased sequence length and sample number, differences in predictors between these viruses would persist. Increases in sample size and sequence length in the smaller and less certain datasets would likely provide clearer and more confident parameter estimates and potentially reveal more subtle effects than could be detected in this study. Additionally, a more complex travel distance modelling strategy that uses a distribution of paths, such as those generated in Circuitscape, may also provide a more realistic measurement of the distance between sites and provide more accurate estimates of these effects (Hall et al. 2021). The increasing accessibility of whole genome sequencing and computational capacity will provide opportunities for future analyses to compare co‐circulating viruses at higher resolution.

In summary, this study showed diversity in the population structures and constraints among six different viruses in a widespread wildlife species, complementing research that suggests levels of host‐pathogen population congruence may vary between systems (Carver and Lunn 2020). Additionally, these results offer insights into possible differences in the biology of these viruses that may be difficult and impractical to measure more directly, such as potential mechanisms of transmission or increased host range. Although viruses can be useful tools for studying hard‐to‐observe host population connectivity, demographics, and behaviour (Biek et al. 2006; Fountain‐Jones et al. 2017), the inconsistencies between bat movement ecology and the connectivity of two common viruses in this study, mastadenovirus and dependoparvovirus, highlight that viruses are not always appropriate tools for reconstructing wildlife movement patterns or the transmission of other pathogens. Viruses used for these studies should be chosen with awareness of their ecology and epidemiology, such as their transmission mode, to establish whether inferences or extrapolations made using them are trustworthy. Similarly, stakeholders in wildlife virus control, whether for conservation, or veterinary or public health interests, should not assume that a new or emerging virus will spread the same way through a host population as do previously studied viruses.

Author Contributions

Conceptualisation: Avery L. Holmes, Laura M. Bergner, Daniel G. Streicker. Formal analysis: Avery L. Holmes. Funding acquisition: Daniel G. Streicker. Investigation: Avery L. Holmes, Alice Broos, Laura M. Bergner. Methodology: Avery L. Holmes, Maxwell J. Farrell, Rita Ribeiro, Daniel G. Streicker. Project administration: Daniel G. Streicker. Resources: Daniel J. Becker, Jorge E. Carrera, Carlos Tello. Supervision: Daniel G. Streicker, Laura M. Bergner, Alice Broos. Visualisation: Avery L. Holmes. Writing – original draft: Avery L. Holmes, Daniel G. Streicker. Writing – review and editing: All authors.

Funding

This work was supported by the Wellcome Trust (217221/Z/19/Z, 218518/Z/19/Z), NSF and BBSRC Ecology and Evolution of Infectious Diseases Program (NSF: DEB 2011069; BBSRC: BB/V003798/1), Leverhulme Trust (PLP‐2020‐362), UK Medical Research Council, and the Edward Mallinckrodt, Jr. Foundation.

Conflicts of Interest

The authors declare no conflicts of interest.

Supporting information

Figure S1: Map of all sample sites and numbers of virus sequences per site, with black dots indicating bat roost sites, bar charts indicating the number of individual virus sequences, with the site name below each chart. Blue and purple dots indicate rabies lineage 1 and 4 sampling sites, respectively. Black outlines show department borders. Table S1: New virus PCRs used for dependoparvovirus and mastadenovirus sequence generation. Table S2: Environmental variables selected for inclusion as potential explanatory variables for virus population structure. Table S3: Variables in the phylogeneticmodels, their sources, and potential interpretations. Figure S2: Biplot for PCA of environmental variables for all sites. The first principal component (Dim1) was comprised mostly of climatic variables, but also some livestock such as horses. PC2 (Dim2) was comprised mostly of human‐associated variables such as livestock, population density, and human footprint. Table S4: Loadings for environmental variables included in the first two principal components (PCs) of a principal component analysis. Variables are ordered by the absolute value of their loadings. Variable references are provided in Table S2. Figure S3: Full phylogenetic tree of mastadenovirus. Bootstrap values are presented for each node. Naming convention follows this format: (BatID_Site). Sequences from the same bat are indicated with additional letters on the ID. Sites are split into departments: AMA = Amazonas, API = Apurimac, AYA = Ayacucho, CAJ = Cajamarca, CUS = Cusco, HUA = Huánuco, LMA = Lima, LR = Loreto. Figure S4: Full phylogenetic tree of betaherpesvirus. Bootstrap values are presented for each node. Naming convention follows this format: (BatID_Site). Sites are split into departments: AMA = Amazonas, API = Apurimac, AYA = Ayacucho, CAJ = Cajamarca, CUS = Cusco, HUA = Huánuco, LMA = Lima, LR = Loreto. Figure S5: Full phylogenetic tree of deltavirus. Bootstrap values are presented for each node. Naming convention follows this format: (BatID_Site). Sequences from the same bat are indicated with additional letters on the ID. Sites are split into departments: API = Apurimac, AYA = Ayacucho, CAJ = Cajamarca, LMA = Lima. Figure S6: Full phylogenetic tree of dependoparvovirus. Bootstrap values are presented for each node. Naming convention follows this format: (BatID_Site). Sequences from the same bat are indicated with additional letters on the ID. Sites are split into departments: API = Apurimac, AYA = Ayacucho, CAJ = Cajamarca, LMA = Lima. Figure S7: Full phylogenetic tree of rabies lineage 1. Bootstrap values are presented for each node. Naming convention follows this format: (SampleID_Department). Figure S8: Full phylogenetic tree of rabies lineage 4. Bootstrap values are presented for each node. Naming convention follows this format: (SampleID_Department). Table S5: Top 10 candidate models in model selection for mastadenovirus presence. DrPV = D. rotundus dependoparvovirus, df = degrees of freedom. Table S6: Top 10 candidate models in model selection for dependoparvovirus presence. DrAdV = D. rotundus mastadenovirus, df = degrees of freedom. Table S7: Variables included in the top models of mastadenovirus and dependoparvovirus presence. CI = confidence interval. VIF = variance inflation factor. Table S8: Parameter estimates for each virus model. The mean estimate value is presented, as well as the upper and lower values of the 80% and 95% highest posterior density intervals. Figure S9: Pairs plot showing relationships between parameter estimates in the model for mastadenovirus. Figure S10: Pairs plot showing relationships between parameter estimates in the model for betaherpesvirus. Figure S11: Pairs plot showing relationships between parameter estimates in the model for deltavirus. Figure S12: Pairs plot showing relationships between parameter estimates in the model for dependoparvovirus. Figure S13: Pairs plot showing relationships between parameter estimates in the model for rabies virus lineage 1. Figure S14: Pairs plot showing relationships between parameter estimates in the model for rabies virus lineage 2.