Papillomavirus in Wildlife

Papillomaviruses (PV) are associated with epithelial malignancies in animals, including cancer in humans. Limited knowledge exists regarding the evolutionary history of non-human PV. We assessed the phylogeography of PV with emphasis in wildlife hosts. We explored the phylogenetic, geographic, and environmental relationships of PV and hosts applying Bayesian inference and spatial analyses of virus and hosts. We found that the available wildlife PV data support previous reports on the higher incidence of fibropapillomatosis over carcinoma in humans and wildlife, being mammals the most common host. We also found geographic bias on the available wildlife papillomavirus (WPV) information towards the Northern Hemisphere, which may have influenced our results to show Europe as the most likely origin of the available WPV lineages. Therefore, we highlight the need for detailed studies on the presence of WPV in regions and species not included in this study (e.g., reptiles from the tropics) to better inform sites of WPV origin, susceptible species, and spillover potential. Future studies of the clinical and subclinical occurrence, distribution, and phylogenetic signatures of WPV may help understand the spread, virulence, and epidemiology of PV in general. From an evolutionary perspective, our overview suggests that WPV are a promising host-pathogen system to untangle questions regarding co-evolution due to its large geographic distribution and occurrence in a wide variety of aquatic and terrestrial wildlife species.


INTRODUCTION
Papillomaviruses (PV) are small, non-enveloped DNA viruses known to produce lesions on the skin (warts) and mucous membranes (condylomas) of various species (Houten et al., 2001;Villiers et al., 2004;Dyne et al., 2018). Cottontail rabbits were the first reported animal hosts of PV (Shope and Hurst, 1933), and rabbit PV was the first DNA virus known to be connected to tumor-growth (Orth et al., 1977). PV hosts are widely varied (Campo, 2002;Wosiacki et al., 2005;Diniz et al., 2009;Herbst et al., 2009;Araldi et al., 2014;García-Pérez et al., 2014), but are mostly represented in mammals and birds (Shah et al., 2010). PV isolates are commonly described as "PV types" (Villiers et al., 2004), some of which are linked to the development of epithelial malignancies and some lineages of cancer in humans, which has highlighted PV's medical and epidemiological importance (Van Doorslaer et al., 2015). Most of the effort to understand the evolutionary history of PV has been focused on human papillomaviruses (HPV) and PVs in domestic animals (Howley, 1986;Villiers, 1989;Houten et al., 2001;Filippis et al., 2002;Diniz et al., 2009;Orlando et al., 2011). Only a few studies are available on wildlife hosts (e.g., Bravo et al., 2010;Gottschling et al., 2011), showing the considerable gaps of knowledge regarding our understanding of PV ecology in the wild. There are over 120 described HPV, from which 30% relate to cancer (Burchell, 2018). Cancer-related HPV are linked to >90% of anal and cervical cancers, ∼70% of vaginal and vulvar cancers, and >60% of penile cancer cases worldwide (Centers for Disease Control and Prevention, 2018). In the United States, ∼41,000 cases of HPV-related carcinoma cases were reported between 2010 and 2015 (Centers for Disease Control and Prevention, 2018). HPV can be categorized as low or high-risk for cancer development. For example, HPV16 and HPV18 lineages have been described as high-risk lineages in humans Ozcagli et al., 2018;Wong and Klapperich, 2018) associated with ∼5% of all the cancer cases worldwide, and ∼3% of all cancers in United States (Jemal et al., 2013;Plummer et al., 2015).
Disease phylogeography explores the association between biogeographic and phylogenetic features of pathogens (e.g., bacteria, virus, protozoa). Phylogeography has been employed in the study of biodiversity (Domínguez-Domínguez and Vázquez-Domínguez, 2009) to understand the interaction between demographic, genealogic, and environmental processes with species-lineages (Avise, 2000;Domínguez-Domínguez and Vázquez-Domínguez, 2009). Phylogeographic approaches are generally used to assess the impact of historical events on the genetic structure of populations, helping to understand speciation events and the configuration and distribution of extant species (Avise, 2000;Vázquez-Domínguez, 2007;Domínguez-Domínguez and Vázquez-Domínguez, 2009). Disease phylogeography has made major contributions to epidemiology via studies of co-evolution between hosts and viruses, in which viruses and hosts reciprocally affect each other's evolution (Webster et al., 2002;Streicker et al., 2010), reconstructions of disease emergence and spread (Mutreja et al., 2011), and assessments of the historical evolution of epidemics (Cohen, 2000;Altizer et al., 2003).
Phylogeographic approaches for public health can be useful to understand the origins, dispersal, and establishment of pathogens through time and populations with reconstructions of transmission. This allows the prediction of future outbreaks via pathogen spillover-i.e., pathogen transmission between different host species (Shah et al., 2010). Geographical barriers can isolate pathogens facilitating the generation of new species (Meinilä et al., 2004) (i.e., vicariance), forcing virus adaptation to available hosts for a long-term relationship stimulating mutual adaptation-i.e., co-evolution (Shah et al., 2010). Ancestral pathogens can evolve into different lineages over time forming new species, strains, serotypes, or other epidemiological classifications. Some of these lineages have varying virulence across host species. For instance, pathogens in original hosts can be primarily avirulent or virulent only under specific conditions such as immunosuppression. Contrarily, severe symptoms can emerge in newer hosts (Morand et al., 2015).
Our understanding of PV phylogeography come from studies focused mostly on HPV lineages and domestic animals, neglecting potential variation in non-human PV (Ho et al., 1993;García-Vallvé et al., 2005;Gottschling et al., 2007;Pimenoff et al., 2016;Van Doorslaer and McBride, 2016). This species-biased research limits our understanding of the evolutionary history of PV and our abilities to predict and control future PV spillover in humans and animals. Limited availability of nonhuman papillomavirus information could be behind the limited research on the phylogeography of non-human PV of wildlife origin. Similarly, the burden and distribution of non-human PVs have not been explored in detail. Therefore, we assessed the phylogeography of PV in wildlife hosts, or wildlife papillomavirus (WPV). This study explored the phylogenetic, host, geographic, and environmental relationships of WPV in a multi-scale ecological assessment. There have been previous studies on the origin and evolution of PV that include non-human species (e.g., Van Doorslaer, 2013;Van Doorslaer et al., 2017;Willemsen and Bravo, 2019). We complement this information by focusing analyses on whole genome sequences from WPV, and environmental, and geographic reconstructions. Our goal was to respond to the questions: Are WPV found under tractable and quantifiable phylogenetic, geographic, and environmental conditions? Can we reconstruct the evolutionary history of the virus? And, can the available data reveal ecological patterns of WPV and its hosts? Exploring these questions can help elaborate new hypotheses of WPV evolution to guide future research in disease ecology and biogeography. A potential new hypothesis may include identification of the origins and paths of dispersal of WPV, which can be assessed using well designed phylogeographic analysis with data not yet available.

Phylogenetic Analyses
We reviewed PV records in animals from scientific literature using the keywords "wildlife" + "papillomavirus" + "genome" in the Web of Science and Scopus databases to collect the most available reports of genetic information worldwide to understand the evolutionary history of the virus. Because our exploration in the Web of Science and Scopus resulted in a limited amount of studies (Web of Science n = 7; Scopus n = 31), we gathered complementary scientific publications in Google Scholar using the same keywords, reviewing studies manually to include only literature matching our inclusion criteria. Our review criteria included PV records with GenBank identification, geographic origin, host species, and reported symptomatology (e.g., asymptomatic, fibropapilloma, carcinoma). Reports without reliable location were removed from the analysis. To better assess the human-wildlife interface (i.e., geographic and temporal overlap between both humans and wildlife, resulting in potential exchange of pathogens (Gortazar et al., 2015), we selected reports with complete genome sequences stored in the GenBank, complemented with PV genomes of human and domestic animals to reconstruct the phylogenetic relationships among wild and domesticated PV lineages.
Sequences were aligned using Mega X (Kumar et al., 2018). Given that the sequences analyzed belonged to the same virus species, we assumed a constant mutation rate and applied the Unweighted Pair Group Method with Arithmetic Mean for the alignment of the sequences, as recommended by the developer (Kumar et al., 2018). Phylogenetic trees were generated by a discrete phylogeography estimation by Bayesian inference through Markov chain Monte Carlo (MCMC), implemented in BEAST v2.5.0 (Bouckaert et al., 2014), through a Hasegawa-Kishino-Yano (HKY) model to account for occurrence of nucleotides at different frequencies, transitions, and transversions occurring at different rates with discrete gamma rate heterogeneity amongst sites (Smith et al., 2009;Sana et al., 2018). Our whole genome matrix was partitioned, following previously described nucleotide positions 1-3702 and 3703-7077 (Gottschling et al., 2007). The statistical uncertainty in parameter values across the sampled tree is given by the 95% highest probability density (HPD) values. All analyses were developed for 200 million generations, sampling every 10,000th generation and removing 10% as chain burn-in following recommendations of the developer (Bouckaert et al., 2014). All the Markov chain Monte Carlo analyses were investigated using Tracer software v1.7  to ensure adequate effective sample sizes (ESS) (>200), which were obtained for all parameters. The final tree was summarized and visualized via Tree Annotator v. 2.3.0 and FigTree 1.4.3, respectively (included in BEAST v2.5.0) (Rambaut and Drummond, 2016;Rambaut, 2017).
To reconstruct WPV phylogeographic radiation, we characterized the geographic distribution of hosts of the reported PV cases considered in this study. We identified the present-day distribution of each host at the continental-level (i.e., Africa, Asia, Europe, North America, South America, Oceania, and Antarctica) due to uncertainties found in the site locations of WPV reports. Host distributions were based on the species' geographic range defined by the International Union for Conservation of Nature (IUCN) (Cosiaux et al., 2016). To assess the most likely geographic origin of the available WPV, we employed Bayesian Binary Markov Chain Monte Carlo (BBM) method analysis to reconstruct the ancestral regions occupied by PV, as previously used in other studies (Durães-Carvalho et al., 2015;Lauron et al., 2015;Forni et al., 2018). BBM was implemented using the software Reconstruction of Ancestral States in Phylogenies, RASP v3.2 (Yu et al., 2015). RASP calculates the probability that a certain parameter representing the origin of a sample was obtained randomly, assessing the quality of the input data. BBM offers a statistical procedure for inferring ancestral states, including geographic distributions at ancestral nodes. BBM assessment uses a full hierarchical Bayesian approach combining the likelihood of the data observed and the probability of events to happen in the order that they did (Ronquist, 2004). We chose to use BBM because it allows null character status information for a portion of input sequences and reduces the uncertainty in the reconstructions of ancestral areas, by reducing the risks of type I error (rejecting the true frequency of monophyletic groups). This attaches a high confidence in large numbers of correct internodes and increases the accuracy of the Bayesian posterior probabilities (Larget and Simon, 1999;Alfaro et al., 2003;Ronquist, 2004).
To account for phylogenetic uncertainty in this analysis, Markov Chain Monte Carlo (MCMC) analysis was conducted simultaneously for one million generations using 10 chains, sampling every 100 generations while the first 25% of generations were discarded as burn-in. Fixed (JC) + G (Jukes-Cantor + Gamma) were used with null root distribution and the temperature was set to 0.1, facilitating the most efficient chain-swapping (Huelsenbeck and Ronquist, 2005). Finally, lineages through time analysis (LTT) was implemented in Tracer  and performed on the previously obtained phylogenetic tree to plot how fast PV lineages have originated over time and assess whether lineage diversification was gradual, or whether it occurred earlier or more recently in time.

Spatial Analyses
To reconstruct the distributions of WPV host species, we revised the geographic ranges of each species for which we found WPV reports according to the IUCN Red List open-access repository (Cosiaux et al., 2016). To assess how many hosts could potentially gather in a specific geographical place, we calculated the host species richness by using spatial analysis in macroecology (SAM) software (Rangel et al., 2010). We created a GIS grid shapefile of the distribution of all hosts species using 50 km spatial scale. We calculated species richness for each taxonomic class to assess species' geographic overlaps. Final maps were created using ArcGIS v10.5 (Environmental Systems Research Institute, 2017).

Environmental Analyses
We used NicheA software to determine the environmental space occupied by WPV and wildlife hosts (Qiao et al., 2016). Specifically, we mapped the distribution of viruses and hosts in environmental dimension. To generate the environmental dimension, we developed a principal component analysis in NicheA based on 19 satellite-derived global bioclimatic variables (2000s decade; 2.5 arcmin resolution) from MERRAclim with information of global temperature and humidity (Vega et al., 2017). The environmental distribution of WPV was estimated using a convex polyhedron calculated around the WPV records and raster of host ranges in terms of the first three principal components (Qiao et al., 2016). We quantified the environmental overlap between WPV and wildlife hosts by taxonomic class using the Jaccard similarity index (Jaccard, 1912). This index is used in ecology to assess the overlap between two organisms in terms of species traits, use of resources, geographic distribution, temporal distribution, and environmental distribution (see Jaccard, 1912;Escobar et al., 2015).

RESULTS
A total of 92 PV records were used to reconstruct the distributions of viruses and host species (Supplementary Table S1). Most samples originated from Europe and North America, presenting notable differences between the numbers of lineages described and their frequency and type of symptomatology. Fibropapilloma was the most frequent form of lesion reported in all continents, and North America and Europe reported all forms of lesions (i.e., carcinoma, fribropapilloma, asymptomatic) ( Figure 1A). Records of asymptomatic hosts were mainly found in Europe. A considerable number of WPV records originated from veterinary reports and captive wildlife species (see Supplementary Table S1). The distribution of WPV hosts clustered in the northern Hemisphere ( Figure 1B), with strong richness of hosts mammals in North America and Europe. Reports of PV in birds occurred in most continents  except Africa and Australia. PVs in bats were reported in Asia and North America, and only North America and Europe reported PVs in rabbits and reptiles ( Figure 1B). Marine and terrestrial mammals were the group of hosts covering most of the geographic ranges used by WPV hosts, although the different classes of hosts presented different geographic ranges (Figure 2). The WPV-infected mammals were reported more often in the northern Hemisphere, with less samples coming from the Neotropic, Africa, and Oceania. Bird hosts gathered mostly in Europe and the coastal areas of the Antarctic. There were only three reports of reptiles, one snake and two marine turtle species, all with broad geographic distributions.
The phylogenetic relationships revealed significant clusters of WPV from hosts sharing evolutionary history (e.g., from the same class). That is, PV lineages found to be phylogenetically related tended to infect hosts from the same clade (Figure 3). However, these associations were not perfect (see Figure 3). The addition of human papillomavirus (HPV16 and HPV18) to the phylogeny revealed a cluster of PV from apes and humans. Ancestral reconstruction and LTT analyses, based on the data available, show that the most likely origin of the common ancestor of the analyzed lineages is Europe, followed by North America and Africa (Figure 4A), from which WPV gradually diverged to the current lineages ( Figure 4B).
Environmental space occupied by reported WPV cases and their wildlife hosts suggested that WPV infections in reptiles are mostly seen in a specific subset of environmental conditions, and that host species in this class occupy territories with environmental conditions beyond those where WPV infections have been recorded (Figure 5). Contrarily, WPV infections in birds and mammals were seen in environmental conditions that overlapped broadly with the space occupied by hosts. This suggests that the environmental range where WPV have been recorded from mammals and birds was similar to the entire environmental range of their host species. These results were supported by the similarity analyses, which revealed that reptiles had a slightly lower similarity observed (Jaccard Index= 0.23) between the environments occupied by the reported WPV and their host as compared with birds and mammals (Jaccard Index = 0.29 and 0.28, respectively) ( Figure 5). Thus, WPV in mammals were found throughout the whole geographic range and environments occupied by their hosts.

DISCUSSION
We assessed eco-phylogenetic patterns of wildlife papillomavirus (WPV). We were able to reconstruct evolutionary, distributional, and ecological signatures of PV lineages using available data to date. Results suggest Europe as the most likely geographic origin of the available lineages. Evolutionary linkages between virus and hosts show linkages between WPV similarity and hosts species infected. The phylogeny of PV did not mirror the host phylogeny, aligned with patterns found in previous studies. The coarse scale reconstruction of WPV distribution identified consistent environmental conditions where WPV can be found, and WPV lineages, geographic locations, and hosts species that deserve future surveillance.
Ancestral reconstruction analyses are widely used to understand the origin and evolution of pathogens (Bargues et al., 2006;Bourhy et al., 2008;Paraskevis et al., 2013;Hayman et al., 2016). Here, available WPV reports allowed us to produce a preliminary overview of the evolutionary history and distribution of papillomaviruses in wildlife species, and to reveal the remarkable gaps of data and surveillance. Our ecological and evolutionary approach explored WPV coevolution and spillover, geographic range of virus and host, and plausible origins that would not be considered in classic epidemiological studies (Galvani, 2003;Faria et al., 2011). A better understanding of the ecology and evolution of WPV may help to anticipate outbreaks in wildlife species and will us help to understand the potential emergence of new pathogenic WPV lineages in humans and livestock (Inhorn and Brown, 1990;Price, 2002;Nasir and Campo, 2008;Trewby et al., 2014).
The ancestral reconstruction revealed a high likelihood of diversification and appearance of new lineages consistently and gradually over time, suggesting that WPV diversification would be linked to host diversification, in agreement with a previous assessment (Shah et al., 2010). This host diversification, along with domestication and the ongoing globalization of PV hosts, could explain the similarity between the homogeneity of genomes collected from different parts of the globe, which poses an analogy to the invasive processes described by Scheele et al. (2019), where they found high genetic diversity in the areas of origin of the studied pathogen but low diversity (i.e. homogenization) on the invaded areas (Scheele et al., 2019). The fact that WPV and their wildlife hosts can be found across broad environmental conditions (Figure 5) suggests that the virus could be infecting wildlife hosts distributed in the southern Hemisphere in areas and species beyond those available to us, at least in the host classes evaluated here (Figures 2, 5). Thus, more efforts are necessary to fill taxonomic and geographic gaps by collecting data from (i) novel species to understand potential WPV ability to infect a wider variety of hosts, (ii) surveillance efforts in neglected systems such as marine, freshwaters, and the Arctic ecosystems, and (iii) sampling neglected locations where known hosts occur, which will facilitate the identification of virus-host codivergence (vicariance) and potential variations in virulence across hosts' populations and geographic regions.
A limited number of wildlife host-species from a broad range of taxa suggests the existence of undiscovered WPV, especially from mammals phylogenetically close to the species included in this study, aligned with reports for other pathogens (Olival et al., 2017). From a clinical point of view, a potential driver of WPV linked to fibropapilloma can be associated to the biased sampling of symptomatic individuals, impeding the assessment of WPV subclinical circulation. It has been previously described that some species are often infected by a different number of PV belonging to a variety of taxonomic groups (García-Vallvé et al., 2005), which, along with the presence of spillover events, can be one of the key factors for WPV transmission success. For example, Nasir and Campo (2008) described in detail how Bovine papillomavirus (BPV) can induce benign tumors in cattle (its main host) but also fibroblastic tumors in equids (Nasir and Campo, 2008), spillover FIGURE 3 | Phylogenetic distribution of WPV lineages in mammals. Dots denote the host species reported according to each WPV lineage. Green: ungulates; purple: rodents; orange: apes: light blue: cetaceans; gray: bats; red: carnivores; dark blue: rabbits (lagomorphs); light pink: humans; yellow: other hosts (FJ379293: hedgehog; AY609301: manatee; GU220391: marsupial). that has lately been corroborated by Trewby et al. (2014) via the analysis of its evolutionary history and cross-species link with equids. However, despite the propensity of sister species to participate in PV spillover, PV phylogenies do not necessarily mirror the phylogenies of their hosts (Chan et al., 1995;García-Vallvé et al., 2005). Thus, a more comprehensive survey of PV lineages in specific wildlife species of broad distributions would increase the understanding on WPV transmission and the likely scenarios of spillover to domesticated animals. As anecdotal note, during the search for WPV records we found clinical and pathological PV reports from diverse wildlife species, however, genetic analyses were generally absent. Routine genetic characterizations of WPV and data sharing is a next challenge for PV research.
The global geographic distribution of WPV (Figures 1, 2), its considerable environmental tolerances (Figure 5), and the enormous range of taxa reported to be infected and symptomatic ( Figure 1B and Supplementary Table S1), support the idea of WPV as an ecological generalist pathogen, making it suitable to study spillover and co-evolution patterns in the animal kingdom.
Other pathogens like rabies and avian influenza have served to understand co-evolutionary forces of viruses and hosts (Streicker et al., 2010(Streicker et al., , 2012Dijk et al., 2015). Thus, PV may be a model system to understand the ecology and evolution of infectious diseases due to its broad range host classes, distributional areas, and environments colonized compared to other viruses. In this sense, WPV can infect a broad range of species including marine and terrestrial mammals, birds, and reptiles across continents and across impressive ecological gradients (Supplementary Table S1). Our exploration may be useful to promote generation of WPV data to better inform veterinary epidemiology, regarding potential emergence of novel WPV in wildlife or livestock and identify scenarios and locations of likely WPV transmission. Enhanced WPV surveillance can also help to elucidate patterns of virus evolution, host adaptation, and virulence. From an ecological perspective, WPV brings an opportunity to study a single pathogen colonizing diverse ecosystems, taxa, and regions, which is urgently needed in disease biogeography to better understand how global change can impact disease transmission.
We found that WPV hosts reported on the available data have a global geographic distribution. Similarly, WPV infections FIGURE 5 | Environmental distribution of the WPV hosts. The first three principal components (X, Y, Z axes) from the original bioclimatic data representing the three-dimensional environmental space (red lines). Distribution of WPV hosts (gray convex polyhedron) and WPV-reports (blue convex polyhedron) are displayed in the environmental space. From left to right: reptiles, birds, and mammals. Gray dots represent single environmental conditions available in the world. observed in mammals and birds were seen in nearly the complete host's range. The capacity of some WPV lineages to overlap almost the entire host distribution suggests the potential of WPV to infect novel species present in environments and geographic areas overlapping the distribution of known hosts. This was supported by finding WPV mammal hosts in aquatic and terrestrial ecosystems, and from birds present in all the continents, including Antarctic regions (Figure 2). We argue that, because multiple WPV can co-infect a single host (Shah et al., 2010), it is possible that sympatric host species could be co-infected with diverse WPV, facilitating the generation of novel papilloma lineages, or that a single WPV lineage could spill over different host taxa and result in a successful establishment of the virus in the new host, as in the case of bovine PV (BPV1) transmission from cattle to horses (Nasir and Campo, 2008;Trewby et al., 2014). This spillover and establishment phenomena has been observed for other zoonotic viruses of wildlife origin in spillover events from wildlife to humans (e.g., Hendra virus and Nipah virus) with establishment of the virus in the naïve species (Plowright and Eby, 2015). These complex and usually unnoticed spillover events could trigger the emergence of novel pathogens through the combination of potential uninfected hosts (humans, wildlife, and livestock), infected host reservoirs, and environmental conditions suitable to maintain the infection (Wu et al., 2018;Johnson et al., 2019). Zoonotic viruses' most common reservoirs are rodents, primates, bats, and domestic animals (Alexander et al., 2018), whose taxonomy and geographical location have been described as the main pillars for these spillover events to occur (Wu et al., 2018), which can also be the main drivers of spillover in WPV.
Our results revealed two forms of potential sampling bias, one expressed as a high number of reports from regions in the northern Hemisphere and another related to symptomatic individuals from studies reporting WPV that produce skin lesions (Villiers, 1989;Villiers et al., 2004;Stevens et al., 2010;Rector and Van Ranst, 2013). We note that our results were restricted to samples available to us and may vary as the number and distribution of available reports increase (Cusimano and Renner, 2010). However, some patterns detected by this study could be a stepping-stone for the assessment of ecological and anthropogenic factors facilitating or limiting WPV spread. Reconstructions of the areas where pathogens could or could not occur have been developed for other infectious diseases, such as Influenza, Filovirus, and Bluetongue viruses, among others (Peterson et al., 2004;Peterson, 2006;Charland et al., 2009;Pioz et al., 2012), facilitating surveillance efforts and pathogen discovery. Importantly, the possibility of the cross-infection between host species (i.e., virus spillover) cannot be ruled out in captive hosts. Thus, our exploration of WPV hosts could guide future surveillance efforts in the wildlife-domestic interface.
An interesting finding is the fact that WPV demonstrates the ability to infect and cause lesions across diverse wildlife species. This highlights the possibility of WPV spreading through naïve populations of the species reported here, which could be of particular concern in the conservation of some species already threatened due to genetic isolation and habitat loss (e.g., apes) (Goossens et al., 2006). This finding may generate new research questions, including whether virulence decreases in the most ancestral host of the virus and whether reducing sampling bias could modify some of the patterns detected in this preliminary overview. Our preliminary maps of WPV distribution can also help identify risk areas for future WPV surveillance, and most importantly, the regions where more surveillance effort is necessary (e.g., tropics). Thus, we argue that a main contribution of this study is the identification of data limitations to identify robust macroecological and phylogeographic signatures of wildlife papillomavirus. Future research should also explore the environmental drivers of transmission, circulation, and virulence of WPV at the local level to determine the most suitable landscape conditions where WPV spillover may occur (Pearman et al., 2008;Peterson, 2014).
Finally, this study represents a preliminary global exploration of the phylogeographic structure of papillomavirus in wildlife, including their geographic distribution to increase our understanding on the recent expansions of virus lineages, the recent appearance of spillover events, and the likely path toward PV origin and evolution. Because it has been suggested that WPV tumors emerge under conditions of stress in the host (Pereira et al., 2003), PV surveillance in free-ranging wildlife may serve as a measure of ecosystem health and as an early warning system of environmental pollution and global change effects over wildlife.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the manuscript/Supplementary Files.

AUTHOR CONTRIBUTIONS
AF-D-D contributed to the design of the project, performed the data collection, the phylogenetic analyses, and wrote the manuscript. MJ performed the spatial and ecological niche analyses and wrote the manuscript. LE contributed to the design of the project, supervised the methodology, and wrote the manuscript.

ACKNOWLEDGMENTS
Authors thank Paige Van de Vuurst and Steven N. Winter for their helpful comments and suggestions. Comments from the reviewers and editor were extremely helpful and greatly improved the content of this manuscript.