DNA virome composition of two sympatric wild felids, bobcat (Lynx rufus) and puma (Puma concolor) in Sonora, Mexico

With viruses often having devastating effects on wildlife population fitness and wild mammals serving as pathogen reservoirs for potentially zoonotic diseases, determining the viral diversity present in wild mammals is both a conservation and One Health priority. Additionally, transmission from more abundant hosts could increase the extinction risk of threatened sympatric species. We leveraged an existing circular DNA enriched metagenomic dataset generated from bobcat (Lynx rufus, n = 9) and puma (Puma concolor, n = 13) scat samples non-invasively collected from Sonora, Mexico, to characterize fecal DNA viromes of each species and determine the extent that viruses are shared between them. Using the metaWRAP pipeline to co-assemble viral genomes for comparative metagenomic analysis, we observed diverse circular DNA viruses in both species, including circoviruses, genomoviruses, and anelloviruses. We found that differences in DNA virome composition were partly attributed to host species, although there was overlap between viruses in bobcats and pumas. Pumas exhibited greater levels of alpha diversity, possibly due to bioaccumulation of pathogens in apex predators. Shared viral taxa may reflect dietary overlap, shared environmental resources, or transmission through host interactions, although we cannot rule out species-specific host-virus coevolution for the taxa detected through co-assembly. However, our detection of integrated feline foamy virus (FFV) suggests Sonoran pumas may interact with domestic cats. Our results contribute to the growing baseline knowledge of wild felid viral diversity. Future research including samples from additional sources (e.g., prey items, tissues) may help to clarify host associations and determine the pathogenicity of detected viruses.


Introduction
Infectious diseases are playing a critical role in wildlife population conservation (Lewis et al., 2017). Through reducing the survival and reproduction of individuals (Woodroffe, 1999;Deem et al., 2001), parasites and infectious diseases can generate trophic cascades (Frainer et al., 2018;Baruzzi et al., 2022), contributing to significant wildlife population declines, affecting multiple species in a community (Pedersen et al., 2007). Furthermore, the transmission of generalist parasites or infectious agents to threatened species can increase the extinction risk of wild animals (Daszak et al., 1999;Woodroffe, 1999;Lafferty and Gerber, 2002;Pedersen et al., 2007). This is evident in carnivores (Woodroffe, 1999;Lafferty and Gerber, 2002;Pedersen et al., 2007), which appear to be particularly susceptible to long term impacts of epizootic diseases at the population-level (Malmberg et al., 2021).
Pathogen surveillance and understanding the viral diversity and potential disease threats associated with more abundant host species may reveal emerging infectious diseases and help prevent future outbreaks and subsequent potential loss of sympatric threatened species. More than just a conservation concern, an increased understanding of the viral diversity in wildlife and the potential for spillover to other species is essential for effectively managing future outbreaks (Olival et al., 2017;Carroll et al., 2018) and is a One Health priority (at the nexus of human, animal and environmental health; Mazzamuto et al., 2022).For example, a study on juvenile and adult red foxes (Vulpes vulpes) in peri-urban areas in Croatia noted the dominant presence of fox picobirnavirus and parvovirus in fecal samples, as well as a novel fox circovirus (Lojkić et al., 2016). With red foxes being the most abundant carnivore in the Northern Hemisphere and the novel fox circovirus being very similar to circoviruses isolated from diseased dogs in USA and Italy, it seems possible to posit that red foxes could serve as wildlife virus reservoirs (Lojkić et al., 2016). Such virome characterizations of carnivores are particularly relevant to advancing One Health priorities, with the order Carnivora being ranked among the top five mammalian orders as a source of zoonotic pathogens (Keesing and Ostfeld, 2021). Monitoring wildlife diseases (and particularly wildlife viromes) using non-invasive fecal sample collection is a nascent field (Pannoni et al., 2022;Mazzamuto et al., 2022;Schilling et al., 2022) and could bridge the gap between passive (e.g., voluntary disease reporting) and active wildlife disease surveillance (e.g., submission of samples from hunted game; Cardoso et al., 2022).
The sociality of a species can also affect the spread of infectious diseases (Sah et al., 2018). In group living or social species, group size was thought to influence disease transmission dynamics (Kappeler et al., 2015;Sah et al., 2018), with larger groups and animals living at higher population densities having higher parasite prevalence and burden (Patterson and Ruckstuhl, 2013;Albery et al., 2020). However, recent research suggests that animals can spatially organize their groups to minimize infections (Albery et al., 2020) and that the interactions within a social group and not only its size drive the spread of infectious diseases (Sah et al., 2018).
Social interactions are not limited to group living, as some relatively solitary species have been shown to exhibit complex social networks with a variety of social partners and interactions (Sah et al., 2018). Thus, despite many felids being solitary, the use of shared environmental resources, as well as occasional conflict or predation within and among felid species, creates opportunities for inter-and intra-specific pathogen transmission.
The Sonoran desert is a unique ecoregion home to four species of wild solitary felids: two being more common, bobcat (Lynx rufus) and puma (Puma concolor), and two listed as endangered, the ocelot (Leopardus pardalis) and jaguar (Panthera onca). Currently, little is known about the exact disease threats and viral diversity associated with these felids. As such, Sonoran desert felids provide both the conservation need and a unique opportunity to assess levels of viral diversity present within and shared between these populations of closely related sympatric host carnivore species. Bobcats and pumas, as the more abundant felids, are easier to sample, and surveys of viral diversity in these species may serve as a proxy for virome characterization or indication of potential viral spillover for the rarer felids.
In this paper, we leveraged a collection of fecal samples of wild bobcat and puma from Sonora, Mexico to determine (1) what DNA viruses are present in wild felids in Sonora and (2) the similarity of fecal DNA viromes between these sympatric species. These data will contribute to the growing field of wildlife viromics and to our understanding of the viral diversity present in wild mammalian species.

Sample collection and processing
Bobcat and puma scat samples were collected from Sonora, Mexico, between 2012 and 2014 ( Figure 1). Scats possibly of felid origin, were collected only if determined to be fresh, based on color, moisture, smell, and texture. Host DNA was extracted using Qiagen's DNeasy Blood and Tissue Kit, and species identification was performed through Sanger sequencing of a region of the mitochondrial cytochrome B gene (Verma and Singh, 2002;Naidu et al., 2011;Cassaigne et al., 2016), as previously described for these samples (Payne et al., 2020). Thirteen puma and nine bobcat scat samples were randomly selected for DNA virome analysis (Payne et al., 2020). Separate viral DNA extraction from scat cross-sections, circular viral DNA amplification, and library preparation followed the protocol in Payne et al. (2020). Sequencing libraries were generated using the TruSeq Nano DNA Sample Preparation kit and sequenced on an Illumina HiSeq 4,000 (2 × 100 bp) at Macrogen Inc. (Korea) in 2018.

Bioinformatics and analyses
The metaWRAP pipeline v. 1.3.2 (Uritskiy et al., 2018) was used to process raw sequencing reads for comparative metagenomic analysis. Within the metaWRAP read_qc module, reads were trimmed using default parameters with Trim Galore v. 0.5.0 (Krueger, 2022) as a wrapper around Cutadapt v. 1.18 (Martin, 2011), human contamination was removed with BMTagger v. 3.101 (Rotmistrovsky and Agarwala, 2011) using the hg38 human genome assembly (GCA_000001405.15), and read quality was assessed with FastQC v. 0.11.8 (Andrews, 2010). Untrimmed reads with human contamination removed were deposited in the Sequence Read Archive. Reads from all samples were co-assembled using metaSPAdes v. 3.14.1 (Bankevich et al., 2012;Nurk et al., 2017) with default parameters outside of metaWRAP, and the scaffolds were then used in the metaWRAP assembly module, allowing for assembly of unused reads with MEGAHIT v. 1.1.3 (Li et al., 2015). We elected to perform one co-assembly across all samples to achieve Frontiers in Ecology and Evolution 03 frontiersin.org better detection of rare taxa and allow for direct comparison of bobcat and puma virome composition. The metaWRAP Kraken 2 module was used to run Kraken 2 (Wood et al., 2019) and generate Krona plots (Ondov et al., 2011) to assess taxonomic composition of contigs in the final assembly, reads from individual samples (subset to 1,000,000 reads), and all pooled reads for each host species, using the premade Kraken 2 viral database (v. 9/8/2022 available at https://benlangmead. github.io/aws-indexes/k2). Contig abundances (in genome copies per million reads) were estimated with Salmon v. 0.13.1 (Patro et al., 2017) using the metaWRAP quant_bins module, and contig taxonomy was assigned using a megablast search against the NCBI nucleotide database (available at https://ftp.ncbi.nlm.nih.gov/blast/db/, downloaded 10/20/2022) using blast v. 2.13.0 (Altschul et al., 1990) outside of metaWRAP for v5 database compatibility, with the output further processed by the metaWRAP classify_bins module for pruning blast hits and assigning taxonomy with Taxator-tk (Dröge et al., 2015).
CheckV v. 1.0.1 (Nayfach et al., 2020) was used to assess quality and completeness of viral contigs. Contigs assigned as viral (and not designated as phages) by classify_bins and which were determined to be complete or high-quality (>90% complete) viral genomes with CheckV were retained for community analyses in R. All downstream analyses were repeated with a second set of viral contigs (>66.7% genome completeness), referred to as our "lower-completeness" set (resulting figures in Supplementary Figures S31-S37). The CheckV quality summary, final taxnomic assignment, top blast hit, and abundance per sample of each contig in the high and lower completeness set can be found in Supplementary Table 1. The R package vegan v. 2.5-7 (Oksanen et al., 2020) was used to conduct viral community ecology analyses on the contigs representing high-quality viral genomes. Alpha diversity metrics (species richness, Simpson Diversity Index, and Shannon Diversity Index) were calculated for each sample, and Wilcoxon rank sum tests were used to determine if significant differences in alpha diversity exist between bobcats and pumas. Beta diversity metrics were calculated among all pairs of samples, both considering contig abundances (Bray-Curtis Dissimilarity Index; abundances were in genome copies per million reads, as estimated by quant_bins) and considering contig presence/absence (Jaccard distance). Kruskal-Wallis rank sum tests were used to determine if beta diversity differed significantly between host species pairs, and Dunn's test was used post hoc to determine which comparisons differed significantly using the FSA R package v. 0.9.3 (Ogle et al., 2022).
Vegan was also used to conduct ordination analyses non-metric multidimensional scaling (NMDS) and principal coordinates analysis (PCoA) with both Bray-Curtis and Jaccard distance matrices, and visualizations were generated using ggplot2 (Wickham, 2016) and ggord (Beck, 2022). To further assess differences in virome composition due to host species, permutational multivariate analyses of variance (PERMANOVA) and analyses of similarity (ANOSIM) were conducted Frontiers in Ecology and Evolution 04 frontiersin.org using Bray-Curtis and Jaccard distance matrices with vegan's adonis2 function and anosim function, respectively. To assess correlations between geographic Euclidean distance and beta diversity, separate Mantel tests were performed using Bray-Curtis and Jaccard distance matrices.

MetaWRAP
We obtained 4,044 contigs of 1,000 to 164,510 nts in length after co-assembly (prior to identification of viral sequences), and Krona plots generated by the Kraken 2 module show the viral taxa represented in the assembly (Supplementary Figure  S1), bobcat reads (Supplementary Figure S2), puma reads (Supplementary Figure S3), and individual sample reads (Supplementary Figures S4-S25). Of the portion of the final assembly matching the Kraken 2 viral database (with contig taxonomy weighted by contig length and coverage), 48% represented viruses within Monodnaviria, with 28% identified as belonging to circoviruses ( Figure 2). Viruses within the family Anelloviridae comprised 16% of the viral portion of the assembly, and phages within the class Caudoviricetes comprised 33%, reflecting viruses likely associated with enteric bacteria of the felids. When reads from each sample of the same host species were pooled (Figure 2), 96% of bobcat viral reads matched to those of the viruses in the family Circoviridae, while puma viral reads largely represented viruses in the families Genomoviridae (51%), Retroviridae (felispumavirus, 14%), Anelloviridae (8%), and class Caudoviricetes (25%). The computational analysis identified 38 complete genomes and 16 additional high-quality viral contigs (using CheckV and classify_bins taxonomy results, not including phages), which were used for downstream analyses. We included an additional 34 medium-quality contigs in our "lower-completeness" set.

Alpha and beta diversity
Species richness, Shannon Diversity Index, and Simpson Diversity Index were the alpha diversity metrics calculated for all samples.
We observed a wider range of species richness values for pumas ( Figure 3A), and pumas had higher median values for each of these metrics (Supplementary Figures S26, S27), although differences in alpha diversity metrics among bobcats and pumas were not significant (richness: p = 0.1677; Shannon: p = 0.1264; Simpson: p = 0.1264). However, using our "lower-completeness" set of viral contigs, we found richness differed significantly between pumas and bobcats (p < 0.05; Figure 3C). Median beta diversity values were greatest between pumas and bobcats and lowest among pumas ( Figure 3B and Supplementary Figure S28), and both Bray-Curtis and Jaccard distances were significantly different among different host species pairings (p < 0.01 for both distances), with significant differences among puma-puma pairings and other host species pairings (p-adj < 0.05 and p-adj < 0.01 with Bray-Curtis and Jaccard distances, respectively, for puma-bobcat pairings and p-adj < 0.05 with both distances for bobcat-bobcat pairings). However, using our "lowercompleteness" set of contigs and Jaccard distance ( Figure 3D), significant differences among host species pairs were explained by puma-bobcat pairings having significantly higher beta diversity than puma-puma (p-adj < 0.01) and bobcat-bobcat pairings (p-adj < 0.05).

Effect of host species and geographic distance
The NMDS and PCoA plots reveal an extensive overlap between bobcat and puma viral communities (Figure 4 and Supplementary Figures S29, S30), although pumas and bobcats with the highest richness levels tended to cluster separately in the NMDS based on Bray-Curtis distances (stress = 0.155), and both PCoA and NMDS (stress = 0.211) based on Jaccard distances revealed puma samples outlying the region of overlap between the two species clusters. The PERMANOVAs revealed a significant effect of host species on community composition using Jaccard distances (p < 0.05, R 2 = 0.09418), but not Bray-Curtis distances (p = 0.497, R 2 = 0.04395). ANOSIM revealed significant differences between host species using both Jaccard (p < 0.05, R = 0.1901) and Bray-Curtis distances (p < 0.05, R = 0.1441). Differences between host species remained significant using the more complex "lower-completeness" set using Jaccard distances (PERMANOVA: p < 0.01, R 2 = 0.12081; ANOSIM: p < 0.01, R = 0.2639) and ANOSIM with Bray-Curtis distances (p < 0.05, R = 0.1513; PERMANOVA: p = 0.649, R 2 = 0.0398). Mantel tests did not reveal significant correlations between geographic distance and beta diversity for either Bray-Curtis (p = 0.813, r = −0.1011) or Jaccard distances (p = 0.165, r = 0.1191).

Discussion
Determining the viral diversity present in wildlife is essential for the management and control of emerging infectious diseases. Owing to the large potential for zoonotic spillover, characterizing mammalian viromes is vital to achieving One Health priorities. In non-invasively collected scat samples from wild pumas and bobcats, we observed diverse circular DNA viruses, including circoviruses, genomoviruses, and anelloviruses. Given that rolling circle amplification (RCA) was performed prior to sequencing, we expected to find a high proportion of circular DNA viruses present (although non-circular DNA is not Bar plots showing taxonomy of reads matching viral database from bobcat and puma samples, as well as taxonomy of contigs matching viral database. Contigs are from the final assembly generated by metaSPAdes and MEGAHIT, and metaWRAP weights taxonomy results of contigs based on length and coverage. Frontiers in Ecology and Evolution 05 frontiersin.org excluded from sequencing). Specifically, one of the complete viral genomes, with particularly high prevalence in two bobcat samples, was identified as Sonfela circovirus 1, the genome of which was originally identified in these two samples (Payne et al., 2020). Additionally, the high proportion of contigs with homology to anelloviruses was also expected, as the diversity of anellovirus genomes isolated from these samples has been previously characterized (Kraberger et al., 2021). Interestingly, we documented the presence of feline foamy virus (FFV) within our dataset as a 4.4 kb contig primarily derived from one puma sample (72% of the reads). This suggests that this is an integrated FFV and was detected through carry-over of host DNA. FFV, a contact-dependent, multi-host adapted retrovirus, is known to cause lifelong infection in both domestic and wild felid species (Linial, 2000;Dannemiller et al., 2020), including pumas and domestic cats (Felis catus; Kechejian et al., 2019;Kraberger et al., 2020). Most studies document prevalence in domestic cats, with comparatively few detecting the presence of FFV in wild felids .
Pumas have shown a high prevalence of FFV and a high frequency of intraspecies transmission in other studies (Kechejian et al., 2019;Dannemiller et al., 2020;Kraberger et al., 2020). Additionally, frequent cross-species spillover of FFV has been documented from domestic cats to pumas due to depredation events , and the presence of the virus in a Sonoran puma may indicate that interactions between wild and domestic felids have occurred. However, it is also possible that this was a result of social spillover from another puma. While FFV is generally considered apathogenic, clinically silent infection has been associated with histopathological changes in domestic cats (German et al., 2008;Ledesma-Feliciano et al., 2019), and further research is needed to clarify implications for feline health.
Our analyses also indicate extensive overlap between bobcat and puma DNA viral communities. The broader diversity of viruses observed in pumas may result from exposure to a wider variety of prey species. Pumas have been observed to predate a broad range of taxa, including ungulates, mesocarnivores, and small mammals (Cassaigne Violin plots showing alpha and beta diversity. (A) Species richness of bobcats and pumas, using contigs representing high-quality or complete viral genomes. (B) Jaccard distances between bobcats, between pumas and bobcats, and between pumas, using contigs representing high-quality or complete viral genomes. Mean Jaccard distance differed significantly between puma-puma pairings and other groups (p-adj < 0.01 for puma-bobcat pairings and p-adj < 0.05 for bobcat-bobcat pairings). (C) Species richness of bobcats and pumas, using viral contigs in the "lower-completeness" set. Mean species richness differed significantly between bobcats and pumas (p < 0.05). (D) Jaccard distance between bobcats, between pumas and bobcats, and between pumas, using viral contigs in the "lower-completeness" set. Mean Jaccard distance differed significantly between puma and bobcat pairings and other groups (p-adj < 0.01 for puma-puma pairings and p-adj < 0.05 for bobcat-bobcat pairings).
Frontiers in Ecology and Evolution 06 frontiersin.org Meyer et al., 2020), whereas bobcats are known to primarily specialize on rodents and lagomorphs (Hass, 2009;López-Vidal et al., 2014;Meyer et al., 2020). Furthermore, apex predators are known to experience greater bioaccumulation of viruses (Malmberg et al., 2021). Pumas are also known to occasionally prey upon smaller felids such as bobcats (Hass, 2009;Prude and Cain III, 2021), and previous studies have documented pathogen transmission from bobcat to puma, putatively through competitive contact or depredation (Franklin et al., 2007;Lee et al., 2017;Malmberg et al., 2021). While such interactions may facilitate viral transmission to Sonoran pumas, the presence of shared viral taxa in bobcats and pumas does not necessarily indicate cross-species transmission. As both pumas and bobcats are known to prey on small mammals, and bobcats have been known to predate deer on occasion (Leopold and Krausman, 1986;McKinney and Smith, 2007), shared viral taxa may instead reflect dietary overlap or shared environmental resources, such as water sources. For example, we identified complete viral genomes matching the rodent anelloviruses Neotofec virus NeonRodL2_5 and Neotofec virus NeonRodL2_6 in bobcats and a partial genome matching Dipodfec virus NeonRodF1_131 (within the phylum Cressdnaviricota) in pumas. These viruses were first isolated from white-throated woodrats (Neotoma albigula) and Merriam's kangaroo rat (Dipodomys merriami), respectively, which could be suitable prey items for both felid species (Meyer et al., 2020). Alternatively, shared taxa that infect these felids may be co-evolved within each species. However, we were unable to determine the strains present within each sample (and species) since contigs were generated by co-assembly. Future research including samples from other sources (e.g., prey items, tissue samples) might help to clarify such host associations. We found that geographic distance among scat samples did not have a significant effect on DNA virome composition at this spatial scale. Our results support previous findings of low levels of spatial autocorrelation of pathogen exposure in pumas and bobcats in Florida, Colorado, and California (Gilbertson et al., 2016). Instead, DNA virome composition appears to be shaped by a combination of host species dependent and independent factors, with extensive virome composition overlap observed between host species using ordination analyses, while PERMANOVA and ANOSIM revealed small yet significant effects of host species. Although pseudoreplication is not suspected, the possibility of some individuals being represented by more than one sample may contribute to the observed effect of host species on DNA virome composition.
Despite providing insight into the possible interactions between host and viral communities, further research is needed to clarify the implications of these results for Sonoran felid health. Fecal virome characterization of non-invasively collected scat samples includes many novel and known viruses derived from the scat depositors as well as prey species and environmental contacts, so the host associations and pathogenicity of each virus in the metagenome is unknown. Of the major viral taxa identified here, viruses in the family Circoviridae may be of most interest in terms of bobcat and puma health. Circoviruses are known as the smallest animal pathogens that replicate autonomously (Fisher et al., 2020). They are found in a number of species [freshwater fish, birds, bats, chimpanzees, minks, elk, and humans (Rosario et al., 2017;Fisher et al., 2020)], although their presence is often subclinical (Fisher et al., 2020). In some birds, circoviruses are considered potentially immunosuppressive (Todd, 2000), suggesting that concurrent co-infections could increase the symptoms and severity of disease (Fisher et al., 2020). Some circoviruses are known to cause clinical disease such as hemorrhagic gastroenteritis in dogs (Anderson et al., 2017;Kotsias et al., 2019), often fatal postweaning multisystemic wasting syndrome in pigs (Chae, 2005;Segalés et al., 2005) and virus (BFDV) in birds [beak and feather disease virus (Todd, 2000)]. The identification of potentially disease-causing circoviruses in bobcats and pumas is of concern for both wild felid population health and conservation. With their A B

FIGURE 4
NMDS plots generated using (A) Bray-Curtis dissimilarity (stress = 0.155) and (B) Jaccard distance (stress = 0.211), using contigs representing high-quality or complete viral genomes. Each point represents the viral community composition within a specific sample. Points are colored by host species, and point size is proportional to species richness of each sample. Ellipses corresponding to the two host species groups are shown at the 95% confidence level. Axes (MDS1 and MDS2) correspond to the two axes of variation.
Frontiers in Ecology and Evolution 07 frontiersin.org propensity for displaying tissue tropism (Todd, 2000), transmission of these circoviruses from more abundant/common species to threatened wild felids could result in catastrophic population declines, particularly if these felids are already immunocompromised from co-infections. Furthermore, the high abundance of Sonfela circovirus 1 reads in two bobcat samples may suggest an active infection. However, with these viruses being identified from non-invasively collected scat samples, further study is needed to clarify host associations and consequences for feline health. Although important feline pathogens present in the scats may have been missed through analysis of the circular viral DNA enriched dataset, these results contribute to the documentation of viral diversity in wild felids. Future studies coupling the characterization of broader virome composition and disease dynamics across sympatric populations of wild mammals could help with the identification of viral threats to wildlife, as well as potentially to humans and domestic animals.

Data availability statement
The datasets presented in this study can be found in online repositories. The name of the repository and accession number can be found at: NCBI; PRJNA922235.