Nucleic acid extraction, library preparation, and sequencing

To extract total DNA and RNA from swab samples, we adhered to a previously outlined protocol [30]. Pooled swab samples underwent homogenization and subsequent filtration through a 0.45 µm polyvinylidene difluoride filter (Millipore, Germany). The filtrate was then subjected to centrifugation at 150,000×g for 3 h. The resulting precipitate underwent treatment with a cocktail of DNase and RNase enzymes to eliminate free nucleic acids. Simultaneously, the QIAamp MinElute Virus Spin Kit (Qiagen, Germany) was employed to isolate the viral nucleic acid.

For the synthesis of first-strand cDNA, the SuperScript™ IV First-Strand Synthesis System (Invitrogen, USA) was utilized along with the reverse primer K-8N, a random primer with a 5′-anchor. The conversion of single-stranded cDNA into double-strand was accomplished using Klenow fragment (NEB, USA) and primer K-8N. The sequence-independent single-primer amplification (SIA) method involved the use of primer K, and magnetic beads (Beckman Coulter, USA) and facilitated the purification of amplification products ranging from 300 to 2000 bp. Sequencing libraries for swab samples were constructed using Nextera DNA Flex (Illumina, USA). Paired-end (2×150 bp) sequencing of the library was performed on the Illumina HiSeq X Ten platform.

Viral contigs assembly and annotation

The raw paired-end sequence reads underwent initial quality control procedures, including the removal of adaptor sequences, primer K sequences, and low-quality reads. The resulting clean reads were then de novo assembled into contigs using MEGAHIT [41]. The assembled contigs were compared using DIAMOND [42] against the nonredundant protein database (nr, download in June 2023) obtained from NCBI (https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz). The E-value cut-off was set at 1E-5 to maintain high sensitivity at a low false-positive rate. The daa file obtained from DIAMOND was introduced into MEGAN [43] to visualize the taxonomy tree of species, and contigs were extracted separately for viruses, each viral family, and genus.

After extracting the sequences, particular attention was given to the conserved genomic regions of each viral group. Notable targets for this classification included the RNA-dependent RNA polymerase (RdRp) gene for RNA viruses, the reverse transcriptase (RT) gene for retroviruses, and the DNA polymerase (DPOL) gene for DNA viruses. Although certain crucial viral entities may not be readily uncovered through broad-spectrum scanning methodologies, the implementation of specific local alignment techniques can serve to augment the depth and precision of viral information retrieval. The ORFfinder tool (https://www.ncbi.nlm.nih.gov/orffinder/) was employed to identify Open Reading Frames (ORFs) in viral sequences. By integrating outputs from both ORFfinder and local BLAST alignments facilitated by the diamond software, the virus sequences procured in this investigation were annotated in accordance with reference genomes. However, genomic RNAs translated through ribosomal frameshifting mechanisms—such as certain virus genera of Astroviridae, Retroviridae, and various other viral families—manual processing were conducted to ensure accurate annotation. Determination of viral hosts was carried out based on their primary hosts using the ViralZone web resource (www.expasy.org/viralzone/). This study is primarily focused on vertebrate viruses; however, the hosts of certain viruses such as the families Circoviridae and Picobirnaviridae remain uncertain. Consequently, viruses with unidentified hosts were only briefly investigated alongside insect and plant viruses. Given that sequences affiliated with coronaviruses and paramyxoviruses have been previously curated and disseminated, this study refrains from offering exhaustive descriptions and analyses of redundant data.

Demonstrating the viral diversity harbored by bat populations in China

Viral diversity relies on conserved genomic features as previously elucidated, and the clustering analysis for viral contigs was performed using CD-hit version 4.7 [44]. The vertebrate virus dataset procured during our research endeavor was combined with publicly accessible bat virus data, with the aim of providing a comprehensive overview of bat-borne viruses prevalent in China. Subsequent to the combination, we partitioned the compiled dataset according to the spatial delineation of provinces, municipalities, and autonomous regions, as well as taxonomic classifications pertaining to the genus of the bat hosts. Following this segmentation, a clustering analysis was conducted. The resulting dataset was analyzed using TBtools-II software to create a heatmap for visualization purposes, as well as to conduct principal component analysis and generate an upset plot [40]. PCA was employed to examine the distribution of virus families across geographic location provinces.

To delve into the diversity of viruses carried by different genera of bats, this study constructed a relational network graph between bat genera and viral families using Gephi software [45]. During the network construction, we employed the Force Atlas layout algorithm to optimize the graphical presentation, enhancing the visual harmony of the network by adjusting the repulsion strength parameter.

Diversity of vertebrate viruses discovered in bats

In the analysis of 378 sequencing libraries, a total of 846 vertebrate virus strains were extracted. By cross-referencing the information with the ICTV report (Table S4), referring to pertinent literature, and computing the p-distance of genes conducive for taxonomic classification, we pinpointed 120 viruses within the vertebrate virus dataset that exhibit characteristics indicative of potential new viral species. Among these newly discovered vertebrate viruses, 61.67% (74/120) were DNA viruses, while 38.33% (46/120) were RNA viruses (Table S2). This subset of novel viruses included 294 strains in total, spanning 16 viral families. In terms of quantity, most of the newly discovered viruses originated from the family Vespertilionidae, followed by the families Hipposideridae and Rhinolophidae. On a per-sample basis, the highest proportion novel viruses were identified in the host family Emballonuridae, with 0.152 novel viruses per sample, followed by Pteropodidae with 0.035 novel viruses per sample, and Vespertilionidae with 0.022 novel viruses per sample (Fig. 1a and Table S2).

Within DNA viruses, adenoviruses exhibited high diversity and a widespread distribution across bat populations (Fig. 2a, Table S18, and S19). Except for two viral families within the order Reovirales, our study’s outcomes closely align with data from other studies within the merged database, indicating a high level of consistency. PCA was performed separately for the sampling regions and the identified virus species in this study. The variables considered were the diversity values of viruses within their respective virus families, in relation to their geographical location and host species. Distinct clustering was observed among the sampling regions, with Guangdong, Guangxi, Zhejiang, and Yunnan exhibiting clear separation. Similarly, among virus families, Picornaviridae, Papillomaviridae, and subfamily Parvovirinae were notably segregated from other families, as illustrated in Fig. 2b and c. Subsequently, statistical analyses were conducted on the host information for each cluster of viral strains identified in the study. The results revealed the number of clusters shared between different bat families and their respective viral families or subfamilies (Fig. 2d, Table S20). Notably, the genera Rhinolophus and Hipposideros, which were previously classified within the same family, harbored the highest number of viruses, comprising 10 viral clusters. These clusters not only included the widely studied RNA viruses but also small, non-enveloped DNA viruses, making these species critical for viral sharing. The second most significant host pair was Rhinolophus and Myotis (mouse-eared bats), which shared 5 clusters. This result is consistent with the habitat overlap observed among these bat species.

Open in a separate window

Fig. 2

Vertebrate-associated state of bat virome. a Diversity of viruses within each bat family. Grouping was performed based on the virus genome types and the source datasets. b PCA at the provincial scale for virus diversity. c PCA at the viral taxonomy scale for virus diversity. d UpSet plot displays the relationships between various bat genera and their associated viral clusters. The lower section lists each bat genus along with the corresponding viral sets. The bar chart above indicates the number of unique viral clusters within each genus. Intersections, representing viral clusters shared among different genera, are illustrated through connecting points and the heights of the bars above, which reflect the number of shared clusters

We examined the relationship between the pairwise p-distances of clustered genes and the actual distances between strains. Most data points followed a specific pattern, showing that combinations with low ratios tended to occur when p-distances were low. This pattern indicates the presence of highly similar viral strains in hosts sampled at distant geographical locations (Fig. S2a). Combinations with a ratio lower than 1 represented the clusters with widespread geographic distribution, with the majority consisting of DNA viruses (18 out of 26 clusters). These clusters included several from the subfamily Parvovirinae and the family Papillomaviridae. RNA viruses (8 out of 26 clusters) were predominantly represented by members of the family Picornaviridae, with only one cluster belonging to Hepeviridae. Combinations with a ratio lower than 1 were marked and displayed on a topographic map, revealing that these clusters are distributed over distances exceeding 1600 km, primarily in plains regions, with relatively fewer occurrences in highland areas (Fig. S2a).

To gain a better understanding of virus spread patterns among these wild mammals, we created a host-virus association network. Network analysis shows that the genera Hipposideros and Rhinolophus occupy central positions, suggesting significant roles in virus transmission. The bat genera Taphozous, Myotis, Tylonycteris, Scotophilus, Rousettus, and Miniopterus exhibit increasingly lower centrality in the network. Picornaviridae, Parvovirinae, and Astroviridae have the highest centrality, followed by the families Caliciviridae, Adenoviridae, Papillomaviridae, and Polyomaviridae (Fig. 3a). Statistical data on the viral loads in bats of varying genera, extrapolated from amino acid levels in conservative viral regions, were collated to generate a comprehensive overview of the spectrum of viruses harbored by bats across diverse regions in China. This data was visually represented to elucidate the viral diversity within different bat genera (Fig. 3a). Notably, bats inhabiting the southeastern and southwestern regions of China demonstrated the highest viral species diversity, in contrast to the relatively low diversity observed in the northern regions and Qinghai province, which may be attributed to factors such as sample size variations or ecological distribution disparities. Intermediate levels of viral diversity are evident in bats from the central and western regions. Analysis of virus carriage by distinct bat genera revealed that members of the genus Rhinolophus exhibit the most extensive viral diversity across most regions in China, with the genera Myotis and Pipistrellus following suit. According to the data provided, the genera in our sample have the representation in terms of sample size.

Open in a separate window

Fig. 3

Distribution of viral diversity within bats in China. a Network plot of bat host genera and the viral families identified in samples. The thickness of the connections between host genera and their corresponding viral families represents the diversity of viruses carried by each host, with thicker lines indicating greater viral richness. A color gradient, transitioning from pink to blue, green, and orange, reflects increasing viral diversity. The length of the connections holds no significance in this model. Node sizes vary, with host genus node size proportional to the diversity of viruses within that genus, and viral family node size reflecting the total viral diversity across all associated host genera. Host genus nodes are color-coded to match their corresponding colors in the text, while all viral family nodes are uniformly colored white for clarity. b Geographic plot of viral diversity per province and pie-charts indicating proportion of each host genera that was sampled

Evolutionary relationships of vertebrate viruses discovered in bats

Among the viral families identified in this study, Picornaviridae and Caliciviridae demonstrate exhibited broad distribution across multiple bat genera, while viruses in the family Papillomaviridae were limited to fewer host species (Fig. 4). Interestingly, while previous studies have identified bat-associated viruses across all four genera of the family Flaviviridae, this study only detected members of the genus Pestivirus in our samples.

Open in a separate window

Fig. 4

Evolutionary relationships of viral sequences within 11 vertebrate viral families. The colored dots denote the viral strains and their corresponding host families identified in this study. Phylogenetic trees were estimated to be using the maximum likelihood method based on viral conserved domain proteins or analysis-used proteins. (Picornaviridae, RdRp; Parvovirinae, NS1; Papillomaviridae, L1; Flaviviridae, RdRp; Caliciviridae, RdRp; Adenoviridae, DPOL; Hepeviridae, RdRp; Astroviridae, RdRp; Retroviridae, RT; Polyomaviridae, LTAg; Orthoherpesviridae, DPOL). Each tree was midpoint-rooted, and the scale bar displays behind the figure title

For other viral families, such as Adenoviridae, Astroviridae, and Hepeviridae, bat-borne viruses displayed host-specific clustering, with the bat-associated Hepeviridae viruses belonging to the genus Chirohepevirus. Novel strains were also identified in the family Polyomaviridae, where viruses were found in multiple genera, including a new branch with recently discovered strains. Within the family Orthoherpesviridae, bat-borne viruses are exclusively found in the subfamilies Betaherpesvirinae and Gammaherpesvirinae, with no bat viruses identified in the subfamily Alphaherpesvirinae. In the case of the family Retroviridae, we detected both endogenous and exogenous viruses. For example, the two Betaretroviruses carried by the lesser dawn bat (Eonycteris spelaea) were likely endogenous, as confirmed by their alignment with bat genomes (Table S5), whereas other retroviruses are likely exogenous.

Notably, we identified novel bat poxviruses in Rhinolophus spp. and other small bats. Although there has been one prior report of a bat poxvirus in China, specifically in eastern bent-winged bat (Miniopterus fuliginosus), where a small segment of the A2L gene was identified (GenBank accession no. KJ651959.1), our study assembled and amplified seven poxvirus contigs from a Chinese rufous horseshoe bat (Rhinolophus sinicus), with a total length of 152 kb. The evolutionary analysis of this virus revealed significant divergence from known poxviruses, with only 60% sequence identity to the closest relative, molluscum contagiosum virus (Molluscipoxvirus molluscum) of the genus Molluscipoxvirus, which can infect humans, suggesting it may represent a new species. The evolutionary relationship of the newly discovered poxvirus was investigated through specific analysis of the gB gene (Fig. 5a). Subsequently, the amino acid sequences of 25 conserved genes were concatenated to construct an evolutionary tree for Rhinolophus bat poxvirus and reference sequences (Fig. 5b). Within the subfamily Chordopoxvirinae, the phylogenetic relatedness of Rhinolophus bat poxvirus is notably similar to molluscum contagiosum virus. Conversely, the bat poxvirus detected in the Americas and Australia displays limited genetic similarity with the strain identified in Asia. Analysis of the major antigenic protein gB revealed distinct antigenicity between these two bat-borne poxviruses. The bat poxvirus identified in this study could potentially signify a new species within the family Poxviridae based on its evolutionary relationship and genetic divergence.

Open in a separate window

Fig. 5

Phylogenetic analyses of viruses of concern. a Phylogenetic tree was estimated using the maximum likelihood method based on amino acid sequence of the glycoprotein B (gB) gene of members in Chordopoxvirinae. b Phylogenetic tree was estimated using the maximum likelihood method based on amino acid sequence alignment of the 25 tandem core genes among Chordopoxvirinae members. Solid triangles are used to mark poxviruses carried by bats found in other continents. c Phylogenetic tree was estimated using the maximum likelihood method based on amino acid sequence of the NS5B region of members in Pestivirus

In addition to the rare detection of poxviruses, a nucleotide segment from a great leaf-nosed bat (Hipposideros armiger)-borne rubella-related virus (Fig. S3d) was also extracted. This sequence may represent a new host for the family Matonaviridae after humans, domestic donkey (Equus asinus), cyclops leaf-nosed bat (Doryrhina cyclops), and yellow-necked field mouse (Apodemus flavicollis) [53]. Importantly, both H. armiger and D. cyclops belong to the family Hipposideridae.

Aquatic-associated viruses were also detected (Fig. S3a, b, c). Two fish astrovirus strains, two Salovirus strains (belonging to the family Caliciviridae), and two Secondpapillomavirinae strains were identified in M. ricketti, which has an unusual piscivorous diet [52]. Fish papillomavirus has been rarely detected in China, and the L1 protein sequences of these two viruses shared less than 40% similarity with all known viruses in this group.

We express our gratitude to Dr. Marco Salemi and Dr. Carla Mavian from the University of Florida for their valuable input on the research methodology.