Viral diversity in wild and urban rodents of Yunnan Province, China

ABSTRACT

Rodents represent over 40% of known mammal species and are found in various terrestrial habitats. They are significant reservoirs for zoonotic viruses, including harmful pathogens such as arenaviruses and hantaviruses, yet knowledge of their hosts and distributions is limited. Therefore, characterizing the virome profile in these animals is invaluable for outbreak preparedness, especially in potential hotspots of mammal diversity. This study included 681 organs from 124 rodents and one Chinese tree shrew collected from Yunnan Province, China, during 2020-2021. Metagenomic analysis revealed unique features of mammalian viruses in rodent organs across habitats with varying human disturbances. R. tanezumi in locations with high anthropogenic disturbance exhibited the highest mammal viral diversity, with spleen and lung samples showing the highest diversities for these viruses at the organ level. Mammal viral diversity for both commensal and non-commensal rats was identified to positively correlate with landscape disturbance. Some virus families were associated with particular organs or host species, suggesting tropism for these pathogens. Notably, known and novel viral species that are likely to infect humans were identified. R. tanezumi was identified as a reservoir and carrier for various zoonotic viruses, including porcine bocavirus, hantavirus, cardiovirus, and lyssavirus. These findings highlight the influence of rodent community composition and anthropogenic activities on diverse virome profiles, with R. tanezumi as an important reservoir for zoonotic viruses.

KEYWORDS:

• Yunnan Province
• rodents
• organs
• anthropogenic disturbance
• zoonotic viruses
• viral tropism

Introduction

Zoonotic viral diseases, such as SARS, MERS and COVID-19, originate from natural transmission of pathogens from animals to humans, posing significant threats to public health. Rodents comprise approximately 41% of mammal diversity worldwide, and serve as natural reservoirs for medically important viruses, including hantaviruses and arenaviruses that cause severe haemorrhagic syndromes in humans across Asia, Eastern Europe, and the Americas [1–3].

Biodiversity hotspots, such as Yunnan Province in Southwestern China, can play a crucial role in understanding spillover dynamics as high diversity and disturbance can create favorable conditions for novel pathogen transmission. In Yunnan, where tourism and farming are common, rodents are known carriers of zoonotic pathogens, including hepatitis E virus [4]. Recent surveys have revealed substantial viral diversity in rodents, including diverse coronaviruses and arenaviruses [5,6].

Yunnan shares cultural and geographical connections to Myanmar, Laos, and Vietnam in the mainland Southeast Asia, an area known as a biodiversity hotspot, including bats and rodents. The geographical range is also an epicentre for the emergence of viruses of significant public health threat, including Hantavirus, Nipah virus and Severe acute respiratory syndrome coronavirus, which are frequently driven by wildlife trade and consumption [7]. In Yunnan Province, Xishuangbanna Prefecture remains a major tourist destination due its high number of attractive natural reserves. The rise in tourism and rubber plantations has led to widespread deforestation. Climate change and habitat degradation alter wild animal populations and increase human-animal contact, enhancing spillover risks. Additionally, habitat fragmentation may lead to non-resident rodent species becoming competent hosts for pathogens. Therefore, understanding the factors contributing to pathogen spillover remains crucial.

In this study, we conducted metagenomic analysis on 681 organs from 124 rodents and one Chinese tree shrew in Yunnan Province. We aimed to characterize viral diversity and abundance in wild and urban rodents, and investigate how the diversity of these viruses vary across landscapes with different human influences. Viral families known to infect humans were characterized based on their genomes or partial sequences.

Materials and methods

Study site and animal sampling

The fieldwork was conducted between July 2020 and March 2021 in Mengla County of Xishuangbanna Dai Autonomous Prefecture, Southern Yunnan Province, Southwestern China. Yunnan is part of Mainland Southeast Asia and borders Myanmar, Laos, and Vietnam. Xishuangbanna has a tropical climate and 50% of the area is covered by forest. Rodents were trapped in areas with varying levels of human activities, including tropical rainforests in Xishuangbanna Tropical Botanical Garden (XTBG), abandoned rubber plantation areas, corn farms in Mengyuan Paradise, and Menglun Town. Foldable metal case traps baited with ripe bananas, with or without peanut butter, were used for trapping, following standardized trapping approaches [8]. Captured animals were anesthetized before being euthanized by cervical dislocation. Organ samples, including the brain, intestines, kidney, liver, lung, and spleen, were collected in RNAlater and stored temporarily in dry ice in the field and at −80°C in the laboratory. To prevent cross-contamination, we maintained strict aseptic conditions by disinfecting equipment, changing gloves, and using sterile instruments as well as a dedicated clean area for sample processing. Morphological identification of captured animals was done by a trained fieldworker, consulting a reference manual (Field Guide to the Mammals of South-east Asia) [9].

Sample processing and library preparation

Animal tissues were homogenized in 1X PBS with stainless beads of 2.4 mm in diameter. A volume of 100 μL of the resulting homogenates was used from each sample to construct pools according to organ, species, and sampling location. Nucleic acids from 200 μL of each pool were purified using the RNeasy Mini Kit (QIAGEN) according to manufacturer recommendations. Purified RNA was quantified using Qubit HS RNA (Thermo Fisher Scientific) and diluted to a concentration of 3.3 ng/μL. The RNA was converted to double-stranded cDNA using ProtoScript® II First Strand cDNA Synthesis Kit and NEBNext® Ultra™ II Non-Directional RNA Second Strand Synthesis Module. cDNA libraries were generated using Twist Library Preparation EF Kit 2.0 for ssRNA Virus Detection (Twist Biosciences, San Francisco, CA) and quantified using Qubit dsDNA BR Assay kit (Thermo Fisher Scientific).

Targeted enrichment and next-generation sequencing

Up to 8 indexed libraries were pooled (8-plex) by equal mass. Twist Pan-Viral Panel v1 (Twist Biosciences, San Francisco, CA) kit containing probes for targeting over 1000 viruses known to cause disease in humans, was used for hybrid-capture. Viral families that are enriched (among others) include: Hantaviridae, Peribunyaviridae, Nairoviridae, Orthomyxoviridae, Arenaviridae, Phenuiviridae, Flaviviridae, Togaviridae, Calciviridae, Picornaviridae, Parvoviridae, etc. [10]. Constructed pools were hybridized to probes for 16 h. Post-capture pools were polymerase chain reaction (PCR) amplified for up to 16 cycles, and the final libraries were quantified and quality controlled using Agilent DNA 2100 Kit (Agilent Technologies). Qualified libraries were paired-end sequenced (PE150) to generate 30 million raw reads on the Novaseq 6000 platform (Novogene, China).

Sequence data processing

We utilized the PIpeline for MetaGenomic Analysis of Viruses (PIMGAVir), a viral metagenomic pipeline, on a cluster for quality control, read filtering, assembly, and primary taxonomic results [11]. To exclude unwanted taxa (Eubacteria, Archaeabacteria, and Plantae), we used the NCBI Refseq non-redundant (nr) database. For the data analysis, we implemented two strategies: assembly and clustering. Specifically, we employed both MEGAHIT and SPAdes for assembly to assess their impact on the number and size of viral sequences in our data [12,13]. Additionally, we generated viral operational taxonomic units (vOTUs) using VSEARCH for clustering [14].

Taxonomic assignment, viral read abundance and diversity

Viral contigs and OTUs were retrieved from PIMGAVir cluster for downstream analysis. Contigs were compared against sequences from the NCBI-nr protein database using DIAMOND v2.0.15.153 (e-value cutoff: 1E−4) [15]. Taxon lineages were obtained with Taxonkit v0.14.0 [16]. Reads were mapped back to contigs and vOTUs using Bowtie2 v2.4.4, and mapped read counts were calculated with SAMtools v1.16 [17,18]. The relative read abundance of viruses consisting of the number of reads for a specific taxonomic group per million of mapped reads (RPM) was calculated using Bioinfokit v2.1.0 [19]. Reads mapped to vOTUs were counted per viral family with a 50-read threshold for inclusion. Shannon Diversity Index was calculated based on viral read abundance of different mammal viral families across locations, rodent species, and organs (Supplementary Method 1) [20,21]. Despite the expected reduction or complete removal of some viral taxa due to our application of targeted viral enrichment, no special accommodation was made in our calculation of the resultant viral diversities. While it is conceivable that genomic diversity as measured without probe hybridization and enrichment may be higher, we note that probe-based hybridization mediates a uniform (non-random) reduction in noise in all samples. In this context, it follows that diversity measures from samples with identical treatments can be compared. Raw reads underwent a new SPAdes assembly, and the generated contigs and PCR barcoding sequences were compared against rodent mitochondrial cytochrome b and cytochrome oxidase subunit I sequences from the MIDORI2 reference database (Supplementary Method 1) [22].

Genome annotation and phylogenetic analyses

Viral contigs from PIMGAVir were annotated using Prokka v1.14.5 30. Confirmation of viral annotation was done via BLAST analysis and detecting putative ORFs using NCBI ORF Finder (https://www.ncbi.nlm.nih.gov/orffinder/). Phylogenetic analysis datasets were constructed based on metadata from the International Community for Taxonomy on Viruses (https://ictv.global/vmr) (accessed on 12 April 2023). Selected viral contigs with > =  200 nt (> =  67 aa) alignment length and no stop codon were used for phylogenetic analysis. Multiple sequence alignment was performed using MAFFT version 7 [23]. Phylogenetic analyses were performed with IQ-TREE multicore version 2.0.3 and ModelFinder (-m MF) used for best-fit model determination [24,25]. Tree nodes were supported by ultrafast bootstrap (1000 replicates).

Recombination analyses

Recombination analysis was performed using SimPlot++ v1.3 [26]. Sequence similarity among selected parvoviruses was assessed using SimPlot basic analysis, and visualized as Network-like structures. Bootscanning was performed with a sliding window of 200 nt and 20 steps, using the Kimura (2-parameter) model and Neighbor-Joining approach with 1000 repetitions [27]. Also, intergenic recombination events were searched using a window size of 100 nt and 1000 permutations.

Statistical analyses

The normality of the datasets and statistical tests were conducted in Rstudio using the rstatix and ggplot2 packages [28,29]. Since none of the datasets exhibited a normal distribution, Pairwise Wilcoxon Rank Sum and Kruskal–Wallis chi-squared tests were used to assess viral read abundance across different parameters [30]. To evaluate the impact of human activities and habitat disturbance among sampling locations, we compared total normalized viral reads from each location. Additionally, we compared normalized reads of each viral family to assess significant differences in abundance levels among rodent species and organs. Furthermore, we analyzed the distribution of viral contigs and vOTUs among various host species and organs. A co-occurrence matrix variable was built by calculating the Pearson correlation coefficients between viral families and threshold for significance was set to 0.05 [31]. To show how viral families cluster together at organ level based on their abundance patterns, hierarchical clustering was performed (Supplementary method 1).

Results

Animal sampling

Between 31st July 2020 and 9th March 2021, we trapped 124 rats and 1 Chinese tree shrew (Tupaia belangeri) across various locations (Figure 1(a, b)). A total of 681 organs were collected, including 125 lungs, 125 livers, 102 spleens, 125 intestines, 102 kidneys, and 102 brains. In human settlements of Menglun Town and corn farms of Mengyuan Paradise, the Asian house rat, R. tanezumi was the dominant species (n = 49). The Indochinese forest rat (Rattus andamanensis) was most abundant in forests, including the tropical rainforest in Xishuangbanna Tropical Botanical Garden (XTBG) and the rubber plantation zone 2.5 km away from Menglun (n = 61). Black rats (Rattus rattus) were captured in Menglun (n = 14). Since only one tree shrew was captured inside XTBG, diversity estimation and statistical analyses focused on rodents.

Figure 1. Sampling sites in Yunnan Province and viral metagenomics. Locations where rodents were trapped are shown in in yellow dots (a). The number of organ samples collected according to rodent species and location (b). Mengyuan Paradise hosts the croplands where rodent trapping took place. Average number of suspected viral contigs and N50 scores for each organ across rodent species (c). XTBG, Xishuangbanna Tropical Botanical Garden.

Metagenomic overview

Forty-six organ pools with an average of 14.8 organs per pool were analyzed using targeted enrichment. Primary processing yielded 128G of clean data with an average of 20,439,059 reads per pool. Suspected viral contig distribution showed a significant difference (p-value = 0.04), with R. andamanensis having the highest and R. tanezumi the lowest contig numbers (Figure 1(c)). The distribution of the suspected viral contigs at the organ level showed the highest values for lungs, intestines, and spleens in R. andamanensis, R. rattus, and R. tanezumi, respectively (Figure 1(c)). A better overall assembly quality was achieved with the kidneys and livers of R. andamanensis, followed by the brains, intestines, and livers of R. tanezumi (Figure 1(c)).

When considering the two assembly data, both methods had similar contig length strength (median N50: 819 bp MEGAHIT, 800 bp SPAdes), but SPAdes generated 26 times more contigs (Supplementary Table 1). A significant difference in vOTU distribution among rodent species was observed (p-value = 0.006), with R. andamanensis showing high vOTU abundance in livers and lungs, R. rattus in intestines and spleens, and R. tanezumi in kidneys and brains (Supplementary Figure 1).

Viral abundance

The metagenomic data identified viruses from over 50 viral families including 15 vertebrate-related viral families. Reads mapped to selected viral contigs were counted and normalized. Adenoviridae, Herpesviridae, Iridoviridae, Retroviridae, Parvoviridae, and Poxviridae showed moderate to high read abundance across organs/rodent species (Figure 2(a, b)). Asfaviridae, Flaviviridae, Calciviridae, and Picornaviridae were found with moderate to low abundance depending on host/organ type. Contigs and vOTUs related to arthropods, plants, bacteria, archaea, and fungi were detected with varying abundance (Supplementary Table 2). Manual scanning of suspected viral sequences confirmed several vOTUs and contigs related to Hantaviridae, Rhabdoviridae, Flaviviridae, Orthomyxoviridae, Peribunyaviridae, and Togaviridae (Figure 2(a/b), Supplementary Tables 3, 4 & 5).

Figure 2. Abundance of selected rodent-associated viral contigs in this study. Sum of viral operation taxonomic units for each indicated viral family (a). Viral read abundance following mapping onto SPAdes (orange) and MEGAHIT contigs (blue) (b). Cooccurrence of indicated viral families across the whole dataset (c). The dark red cells in the heatmap represent strong positive correlations between viral families, indicating that these families tend to co-occur together frequently in the dataset. Hierarchical relationships between rodent organs based on the indicated viral family abundance patterns.

The assessment of potential co-occurrences of the selected viruses based on clustering of the dataset showed that pairs of viral families Paramyxoviridae/Alloherpesviridae, Parvoviridae/Calciviridae, Picornaviridae/Calciviridae, Peribunyaviridae/Hantaviridae, and Picornaviridae/ Parvoviridae, tended to co-occur in our dataset (Figure 2(c)). In the context of viral abundance patterns across organs, a distinctive distribution pattern was observed in rodent kidneys (Figure 2(d)). However, comparable abundance patterns were evident among lungs and spleens, as well as brains and intestines (Figure 2(d)).

Viral diversity, association with organ, host species or location

Remarkable mammal viral diversity was found in human-modified landscapes among sampling locations, with R. tanezumi and R. andamanensis showing similar profiles (Figure 3). Mammal viral diversity was found to be positively correlated with the level of habitat disturbance. R. tanezumi in human settlements had the highest viral diversity, while R. andamanensis in rubber plantations exhibited relatively moderate diversity of these viruses (Figure 3(a–c)). R. andamanensis, with habitat range limited to forests and rubber plantations showed an overall lower mammal viral diversity. Interestingly, when considering the two habitats for this rodent species, we found that individuals from rubber plantations exhibited relatively high mammal viral diversity compared to those from the tropical rain forest (Figure 4(a)). Also, R. rattus showed a low mammal viral diversity. Significant differences were found between organs for viral families Calciviridae, Hepeviridae, and Parvoviridae (p-values < 0.05). Parvoviruses were notably present in brains, kidneys, and spleen samples. A high degree of differences in viral read abundance was also observed between host species for Flaviviridae, Iridoviridae, and Picornaviridae (p-values < 0.05).

Figure 3. Viral read abundance according to organ, species, and location. Menglun (a); Mengyuan Paradise (b); Rubber plantations (c); Xishuangbanna Tropical Botanical Garden (d). For each location, the first heatmap (on the left side) shows normalized reads from mapped reads to SPAdes contigs, and the second heatmap displays normalized reads from mapped reads to MEGAHIT contigs. The white spaces in (a) & (d) refer to missing organs for the corresponding animal species.

Figure 4. Overall mammal viral diversity in each rodent species across locations with varying influence of human activities and habitat disturbance (a). Mammal viral diversity of each rodent species at the organ level (b).

Assessing mammal viral diversity by location and organ showed variations in Shannon diversity within and between each of these factors (Figure 4(a, b)). Spleen and lung samples consistently exhibited high diversities of these viruses in all three rodent species. R. tanezumi and R. rattus had similar mammal viral diversity profiles across organs, with lower diversities in brain, kidney, and liver samples of R. rattus. Spleen samples in commensal rodents displayed the greatest mammal viral diversities, followed by lung samples, while brain samples had the lowest diversity for these viruses. R. andamanensis showed relatively lower mammal viral diversities in kidney, intestine, and liver samples.

We further assessed associations between viral read abundance and host organ or species. The comparison of virus abundance between organs revealed significant differences for viral families Calciviridae (p-value = 0.01), Hepeviridae (p-value = 0.001), and Parvoviridae (p-value = 0.01). Parvovirus reads were notably found in brains, kidneys, and spleen samples (Figure 3(b and d)). Considering viral read abundance and host species, a significant difference was found for viral families Flaviviridae (p-value = 0.02), Iridoviridae (p-value = 0.03), and Picornaviridae (p-value = 0.03). Parvoviruses were particularly found in R. tanezumi in Menglun Town and Mengyuan Paradise (p-value = 0.04).

Characterization of viruses

Parvoviridae: An extensive number of parvovirus contigs and OTUs belonging to the subfamily Parvovirinae were identified by both assembly approaches in this study. Contigs related to a new protoparvovirus, tentatively named Banna rat parvovirus 1 (BRPV1, OQ878961), were found in one brain and one kidney pool of R. rattus. The genome of BRPV1 consists of 4030 nt with four ORFs including 2 non-structural proteins (NS1, 467 aa; NS2, 153 aa) and 2 overlapping structural proteins (VP1, 726 aa; VP2, 587 aa). BRPV1 contigs from both assembly methods were found to be generally identical (∼99.9%), and the closest relative is Mouse parvovirus 4b (FJ445512) with a nucleotide identity of 86.8%. Bamboo rat parvovirus and Mouse parvovirus 4b are related to BRPV1 based on NS1 and VP1/VP2, with amino acid identities of 94.5% (AXR9642) and 91.2% (ACJ63478), respectively (Table 1).

Table 1. Pairwise amino acid identities in NS1 and VP1-VP2 between parvoviruses identified in this study (in italics) and indicative parvoviruses.

Parvovirus	Accession no.	Genome / contig size (nt)	NS1		VP1-VP2
			Length (aa)	Identity (%)	Length (aa)	Identity (%)
BRPV1	OQ878961	4030	467	100	726	100
Mouse PV 4b	FJ445512	4794	672	93.0	729	91.2
Bamboo rat PV	MF497825	4754	672	94.5	541	60.5
Rat bufavirus SY-2015	KT716186	4634	644	47.6	721	36.1
Porcine bufavirus	MK279312	4189	579	45.8	679	38.2
Canine PPV	MH645362	4249	638	47.0	704	36.2
Megabat bufavirus 1	LC085675	4765	642	47.4	719	35.7
BPBoV1	OQ878963/2*	458/1586*	149	100	523	100
Porcine bocavirus	MG365887	4798	601	100	601	99.6
Porcine bocavirus	KJ622366	4630	635	97.3	616	92.9
Porcine bocavirus	KY426738	5087	635	96.6	616	92.9
Bocavirus pig/SX/China/2010	HQ223038	4786	636	90.6	623	79.7
Porcine bocavirus 1	HM053693	5153	703	69.1	705	51.5
Porcine bocavirus 3	JF512472	5082	800	65.1	661	38.1

Notes: aa, amino acid; nt, nucleotide; BRPV1, Banna rodent parvovirus 1; PV parvovirus; PPV, protoparvovirus; BPBoV1, Banna porcine bocavirus 1.

*NS1/VP1-VP2.

Notably, several contigs ranging from 209 to 1586 nt related to a porcine bocavirus named Banna porcine bocavirus 1 (BPBoV1) were found in at least one sample of each rodent species. The prevalence of BPBoV1 was higher in R. tanezumi with 4 pools implicated [brains (2), liver (1), kidney (1)], and R. andamanensis [spleen (1), kidney (1)] and R. rattus [brain (1), kidney (1)] each have 2 positive pools. The longest contig (1586 nt, OQ878962) of BPBoV1 was revealed to be the VP1/VP2 ORF. Comparison of the VP1 and VP2 sequence of BPBoV1 to the sequences from the NCBI nr database showed that porcine bocavirus (strain YY.22 isolated from R. norvegicus in China) was the closest hit with nucleotide and amino acid identities of 98.9 and 99.6%, respectively. All contigs of BPBoV1 were further mapped to the genome of YY.22 (MG365887) and 458 nt (149 aa, OQ878963) consisting of a partial NS1 were obtained. The NS1 sequence of BPBoV1 shares 99.56% and 100% of nucleotide and amino acid identity to YY.22 (Table 1).

Phylogenetic analyses placed BRPV1 and BPBoV1 in the genera Protoparvovirus and Parvobocavirus, respectively (Figure 5(a–d)). BRPV1 clustered together with Munite virus of mouse (J02275) and Rat bufavirus SY-2015, based on the whole genome and the NS1 protein, respectively. However, considering the VP1 and VP2 protein, BRPV1 is a root of a branch including Rat bufavirus SY-2015 and other parvoviruses of the genus Amdoparvovirus. Phylogenetic analyses based on VP1 and VP2 positioned BPBoV1 in a separate branch with Porcine bocavirus H18 and Bocavirus pig/SX/China/2010 isolated from pigs in China. Taking into account the proposed species demarcation criteria of <85% amino acid sequence similarity, BRPV1 and BPBoV1 have been identified as Munite virus of mouse and Porcine bocavirus H18 species, respectively.

Figure 5. Maximum likelihood phylogenetic trees of the whole genome (a), NS1 protein (b), and VP1/VP2 (c & d) protein of representative members of the subfamily Parvovirinae. The blue and green colours highlight the genera Bocaparvovirus and Protoparvovirus, respectively. The rodent protoparvovirus (BRPV1, OQ878961) and the procine bocavirus (BPBoV1, OQ878962) identified in this study are shown in orange. The trees were inferred with IQ-TREE multicore version 2.0.3 with GTR + F + R6 (complete genome), LG + F + I + G4 (NS1, BRPV1), Q.pfam + F + I + G4 (VP1, BRPV1), and LG + F + I + G4 (VP1/VP2, BPBoV1) substitution models following alignment with MAFFT version 7 [23, 32]. The phylogenetic tree was supported by Ultrafast bootstrapping with 1000 replicates. The scale bar represents the number of substitutions per site. BRPV1, Banna rat parvovirus 1, BPBoV1, Banna porcine bocavirus 1.

Calciviridae: A few contigs related to a norovirus were recovered in one kidney pool of R. tanezumi (Supplementary Table 5). BLAST analyses supported partial sequences of viral capsid (VP1) and nonstructural protein (NSP) of a novel Murine norovirus named Rodent norovirus 1 (RNV1). RNV1 shares 76.1% (NSP) and 74.4% (VP1) amino acid identity to Murine noroviruses of genotype V (AFV48050, AFV48050) detected in R. norvegicus and Apodemus sylvaticus in Hong Kong and the United Kingdom. Representative members of the family Calciviridae were used to infer the phylogenetic relationship of RNV1 based on partial NSP protein (78 aa). Based on the partial NSP, RNV1 was identified with Murine norovirus 1 at the root of a clade hosting other noroviruses, including Bat norovirus, Norovirus dog/GVI.1/HKU_Ca026F/2007/HKG, and Norwalk virus (Figure 6(a)).

Figure 6. Maximum likelihood phylogenetic tree for a partial nonstructural polyprotein of the family Calciviridae (a) and Picornaviridae (b). The viruses identified in this study are shown in orange/red. The tree was inferred with IQ-TREE multicore version 2.0.3, using the Q.pfam + I + G4 (RNA) and Q.pfam + F + R3 (Banna rat cardiovirus 1) substitution model following alignment with MAFFT version 7 [23,32]. The phylogenetic tree was supported by Ultrafast bootstrapping with 1000 replicates.

Picornaviridae: A new picornavirus, named Banna rat cardiovirus 1 (BRCV1), was recovered from the same kidney sample of R. tanezumi as RNV1. Contigs and vOTUs related to BRCV1 were found (Supplementary Table 5). BRCV1 shares 94.3% and 85.7% of amino acid sequence identity with cardioviruses isolated from a wild rat (MH976711) in China and a human (AWK02670) in Latin America. Phylogenetic analysis placed BRCV1 closer to Saffold virus 1 (cardiovirus D1) (EF165067), which was isolated from a paediatric patient with fever of unknown origin (Figure 6(b)) [32]. BRCV1 also shows a close phylogenetic relationship with other cardioviruses associated with diverse clinical symptoms in rodents (X56019) and humans (AAA43037).

Viruses with low representation

Six other viral families (Hantaviridae, Rhabdoviridae, Flaviviridae, Peribunyaviridae, Orthomyxoviridae, and Togaviridae) with low representation were also noteworthy. It is important to note that the sequences related to these viruses are relatively short (< 400 bp). However, these sequences consist of partial ORFs that were manually confirmed after being confidently detected using DIAMOND BLAST search against the NCBI-nr protein sequence database (Supplementary Tables 4 & 5).

Bunyavirales: Evidence of Bunyavirales (Hantaviridae/Peribunyaviridae) was found, encompassing sequences ranging from 140 to 381 bp. Hantaviruses and orthobunyaviruses are characterized by a tripartite genome, consisting of a long, a medium, and a short segment that encode RNA-dependent RNA polymerase (RdRp), glycoprotein (GP) precursor, and nucleoprotein (NP), respectively. Specifically, we identified three contigs (one RdRp, two NP) and nineteen vOTUs (eight RdRp, five GP, and six NP) associated with Orthohantavirus (Hantaviridae), as well as five vOTUs (four GP, one RdRp) related to Orthobunyavirus (Peribunyaviridae). These viral sequences were identified in R. tanezumi, with hantavirus sequences distributed across various organs (brain, kidney, liver, intestines) compared to orthobunyavirus sequences, which were primarily localized in the brain. Notably, the hantavirus sequences exhibited close genetic relatedness to Hantavirus MAYOV (RdRp, AMP83216), Alder hantavirus (GP, KP013579), and Cao Bang virus detected in Vietnam (NP, EF543524) (Supplementary Table 4). The sequences related to orthobunyavirus displayed substantial amino acid similarity to Douglas virus (GP, QCT81303) and Weldona virus (RdRp, QLA46879), in their respective GP and RdRp genes (Supplementary Table 4).

Rhabdoviridae: Within the family Rhabdoviridae, we found evidence of lyssavirus sequences (Supplementary Table 4). This includes three vOTUs ranging from 149 to 244 bp, all of which were isolated from a single brain pool of R. tanezumi. Through comparative analysis of both nucleotide and amino acid sequences, these rhabdovirus sequences demonstrated a substantial degree of relatedness to the RdRp gene of European bat 1 lyssavirus (CAG9051531) and Taiwan bat lyssavirus 1 detected in the brain of Pipistrellus bats in Taiwan (CAG9051531).

Orthomyxoviridae: We identified a total of 13 viral vOTUs of influenza A virus with a length of 149-309 bp. These sequences corresponded to partial segments of the influenza A virus genome, specifically encoding for the proteins basic 1 (PB1), basic 2 (PB2), and acid (PA). These sequences were detected in four organ samples including brain, kidney, liver, and intestines of R. tanezumi. Notably, the PB1 segment of the Influenza A virus – A/Hong Kong/516/97(H5N1) (AAK49370), the PB2 segment of Influenza A virus (ADC34530), and the PA segment of Influenza A virus (QTK13600) exhibited the closest sequence homology to the Influenza A virus sequences uncovered in this study (Supplementary Table 4).

Flaviviridae: Flaviviruses were identified in two rat species, including R. tanezumi (lung, intestines, brain, and kidney) and R. rattus (spleen). This consists of 5 contigs and 35 vOTUs, with sequence lengths ranging of 149-310 bp. Notably, when compared to sequences in the NCBI database, the best match was found with Langat virus (AAF75260), a member of the Flavivirus genus (Supplementary Table 4).

Togaviridae: Sequences related to the Togaviridae family were also identified in R. tanezumi (brain, lung, intestines). These sequences with 149-275 bp in length, included 2 contigs and 7 vOTUs and displayed a strong similarity in the NSP and capsid protein (CP) regions to alphaviruses. Notably, the top hit in NCBI searches corresponded to Eilat virus (QBG67156.1), an alphavirus primarily associated with insect hosts [33]. The alphavirus sequences also exhibited significant amino acid similarity to the CP of Chikungunya virus (ADO34948), further emphasizing a connection between the identified alphavirus sequences and previously documented alphaviruses (Supplementary Table 4).

Recombination analyses

We explored putative recombination candidates in the BRPV1 genome using sequence network similarity and bootscanning (Figure 7). High sequence similarity was found between BRPV1, Bamboo rat parvovirus, and Mouse parvovirus 4b, suggesting possible intragenic recombination events (Figure 7(a/b)). Notably, two recombination candidates were detected between BRPV1 and Bamboo rat parvovirus (nucleotide positions 1400–1780) and BRPV1 and Mouse parvovirus 4b (nucleotide positions 2680-2980) (Figure 7(c)). Similarly, strong sequence similarity was found between Canine protoparvovirus, Rat bufavirus, and Porcine bufavirus, suggesting potential recombination events among parvoviruses in nature (Figure 7).

Figure 7. Sequence similarity network and recombination analysis of BRPV1. The Simplot results display the variations between BRPV1 and the generated consensus sequence, and similarity network shows the global (grey lines) and local (red lines) similarities between BRPV1 and the indicated parvoviruses (a and b). Bootccanning shows the magnitude of clusters between BRPV1 and other parvoviruses (c). The percentage of permuted trees in the ordinate supports the clustering according to each sliding window. C1 and C2 display nucleotide positions where supported recombination breakpoints in the genome of BRPV1 were noted.

Discussion

Over 60% of emerging infectious diseases come from zoonotic pathogens, with over 70% of these originating from wildlife [34]. Novel or known highly virulent viral diseases pose increasing threats to global health and the economy, causing widespread illness and death. Rodents, being the most diverse mammal group, are widespread and some Muridae species are known reservoirs for notable viruses [4,35]. Yet, many rodent diversity hotspots remain almost unsampled, such as the distribution of such viruses in Yunnan's rodent population. Utilizing high-throughput sequencing technology, we focused on characterizing viral diversity and its association with their rodent hosts and habitat disturbance. Distinct viral traits were identified in rodents from regions with varying human activity and habitat disturbance. Rattus tanezumi in habitations and agricultural fields exhibited high diversity of mammal viruses, suggesting its potential role in carrying zoonotic viruses. The results suggested that rodent virome profiles are influenced by the composition of rodent communities across varying human activities, suggesting further studies on the commensal rat R. tanezumi would be useful for disease preparedness.

Our initial analysis of rodent-associated viruses provides valuable insights into viral ecology and evolution. Despite a limited sample size of 124 individuals, this study offers a crucial glimpse into the viral communities of rodents in Yunnan. Future investigations, incorporating serological or PCR screening, hold the potential to provide more in-depth insights into past infections and viral dynamics. It is worth noting that some viral taxa may have evaded detection either as a result of targeted enrichment for mammalian viruses, or due to other limitations associated with the probes used. To capture a more complete diversity of all viruses present (in addition to mammalian viruses), a non-targeted approach would be more appropriate. Additionally, an extended, year-round sampling would capture seasonal variations more effectively.

Applying different bioinformatics tools to analyze metagenomic data enhances understanding of animal virome. Assembler choice depends mainly on scientific questions, resource availability, and bioinformatics expertise. In this present study, SPAdes showed better viral genome assembly performance, providing larger contigs and higher N50 values across organ samples. Previous investigations also supported this evidence, while MEGAHIT was identified as a resource-saving alternative [36].

We identified co-occurrence patterns among several viral families. Such associations are mainly due shared hosts, evolutionary similarities, and ecological factors like habitats. Unique viral abundance patterns were also observed in different organs, emphasizing the importance of considering the specific environment of each organ. Organs may vary in their specialized receptors, immune responses, and microbiota composition, all of which influencing virus-host interactions as well as the abundance and the co-occurrence of viral families. Understanding these patterns is vital for deciphering viral ecology within hosts.

Our work on viral diversity differs from previous investigations due to distinct study approaches [4,37]. While metagenomics and metatranscriptomics have previously characterized rodent virome at organ level, much remains to be understood regarding the diversity of zoonotic viruses at species and habitat levels [35,38]. This knowledge gap underscores the complexity of host-virus interactions in Yunnan rodent populations, with unique viral profiles across species and organs, with spleen samples consistently exhibiting high mammal viral diversity.

A key question we tried to address in this study is the mammal viral diversity and abundance in relation to rodent species and the degree of anthropogenic disturbance. This is an important question as generalist rodent species, which are most likely to co-exist with humans and animals in close proximity in highly disturbed habitats are frequently more competent hosts than species dependent on intact habitat. Here we found that even within host species from different habitat types such as R. tanezumi, the highest mammal viral diversity was found in the most disturbed (urban) environments. This is consistent with previous studies demonstrating that more competent hosts are often more prevalent in disturbed environments, in which the movement of these animals across a gradient of habitat degradation can potentially introduce pathogens into less disturbed settings [39–42]. Our study revealed that land use change may potentially increase pathogen transmission risk by enhancing the proliferation of commensal rat associated with an increased number of zoonotic viruses [41].

Environmental factors driving the spillover of viruses in rodent populations include, among others: climate change, land use change, land fragmentation, and human expansion into natural habitats [43]. While extensive research on the environmental drivers of zoonotic viral diseases has traditionally focused on climate-related factors, there is a growing consensus that alterations in land use patterns, such as the conversion of natural settings into agricultural, urban, or other human-altered ecosystems, play a significant role in mediating infection risk in reservoir populations and the emergence of spillover events [41,44]. Our results showed even R. andamanensis (a non-commensal rat) individuals found in relatively disturbed areas (rubber plantations) harboured high mammal viral diversity compared to those from low disturbed habitats (tropical rain forest). Moreover, while R. andamanensis was the most sampled rodent species, its low mammal viral diversity was in accordance with its habitat range restricted to landscapes with low degree of human influence. This aligns with previous research findings that have identified areas with significant land use changes in biodiversity hotspots as high risk areas for the emergence of zoonotic viral diseases [45]. For example, several hantavirus strains causing pulmonary syndrome and Brazilian haemorrhagic fever virus (Arenaviridae) were found to be located in fragmented forests and agricultural lands in Brazil [45–48]. In Africa, many biodiversity hotspots were also identified overlapping endemic areas of Lassa virus (Arenaviridae) and mpox virus (Poxviridae) [45,49,50]. Understanding these landscape-level patterns in viral diversity contributes to our knowledge of viral ecology and its response to human influences.

Evidence of similar mammal viral diversity profiles across organs among closely related rodent species (R. tanezumi and R. rattus) was found. This suggests that these two Rattus species have potentially comparable characteristics including the synanthropic behaviour and factors related to susceptibility to viral infection. Moreover, certain viral families exhibited strong associations with specific organs and rodent species. For instance, sequences related to several relevant viral families were predominantly found in the brains of R. tanezumi. This implies that many of these viruses may have tropism for rodent brains. These findings have significant implications for our understanding of zoonotic potential, necessitating further investigation into the distribution of these viruses within wild rodents, including specific organs and tissues.

We identified similar genetic characteristics between the major components of porcine bocaviruses (PBoVs) in our study and newly discovered group 4 PBoVs found in commensal rodents including R. norvegicus and R. tanezumi in China [6]. This discovery confirms the presence of group 4 PBoVs in Southwestern China and suggests the need for further investigation into their potential spread within the country. The identification of R. tanezumi in human settlements, where backyard pig breeding is common, supports the idea that porcine bocaviruses could be transmitted between R. tanezumi and pigs, highlighting the role of these commensal rats as potential carriers of zoonotic viral diseases.

Natural recombination within parvoviruses has been observed in animals, such as various bat species, and mink (Neogale vison) [51,52]. The discovery of a new type of rodent protoparvovirus (BRPV1) in the brains and kidneys of R. rattus showed how these viruses may have evolved and the possibility of genetic recombination events occurring amongst these pathogens. Interestingly, the presence of recombination breakpoints in critical gene regions of BRPV1 implies fascinating instances of genetic exchange within these viruses, providing insights into how parvoviruses have evolved alongside their rodent hosts.

The identification of rodent norovirus and cardiovirus in kidney samples of R. tanezumi points to organ-specific viral preferences for these Picornavirales. Notably, human noroviruses can cause gastroenteritis. The identification of human noroviruses in wild rodents suggest the potential these animals may have for viral spread into the environment [53]. Similarly, the cardiovirus detected in this study shares a strong phylogenetic relationship with viruses known for causing clinical symptoms in humans [54–56]. Further research is imperative to comprehend the evolution and spillover risk of rodent cardioviruses.

The present study extends previous investigations by detecting both hantavirus and lyssavirus in a commensal rat, R. tanezumi. While rodents are natural reservoirs for hantaviruses, the detection of rabies virus in wild rodents is uncommon globally. However, rabies virus has been detected from few wild rodent species in three Chinese provinces [57]. Moreover, rat-associated human rabies cases have also been documented in China [57–59]. The role of commensal rats in Yunnan as potential rabies spreaders to domestic animals, pets, or to humans requires further study.

The identification of animal viruses has been predominantly focused on screening non-invasive samples, including faeces, respiratory tract, and oral cavity from animals in their natural habitats. To assess the potential threat and provide reliable information about the animal source, detecting viral nucleic acid in tissue or blood samples is essential. Using diverse organ types, our study offers insights into mammal viral diversity, abundance, and possible viral tropism in an important animal group.

We provide a comprehensive overview of the mammal viruses with most being human-related viruses associated to rodents across various habitats, shedding light on the diversity, distribution, and potential zoonotic implications of various viruses. These findings emphasize the importance of improving surveillance approaches to better understand the complex interactions between rodents, viruses, and their environment, with potential implications for public health and wildlife conservation.

Ethics approval and consent to participate

The procedure used for sampling rodents was approved by the ethics committee of XTBG, Chinese Academy of Sciences. Permit for trapping and sampling rodents in natural reserves outside XTBG was provided by the local authorities of Xishuangbanna Natural Reserve (Approval number: 384).

Author's contributions

Conceptualization: GW and NB. Design: GW, NB, and YK. Sample collection: AT, YK, RL, YC, and ACH. Sample processing, library preparation and NGS: AT, NB and YK. Bioinformatics and data analysis. YK, NB, EM, GW, AT and JL. First draft of the manuscript: YK. Manuscript review and editing: GW, NB, ACH, and AT. All authors have reviewed the manuscript and agreed to the submitted version for publication.

Acknowledgments

We would like to express our gratitude to the Landscape Ecology Group at Xishuangbanna Tropical Botanical Garden for hosting the research team during the fieldwork portion for this project.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Sequencing data generated during study was submitted to China National Microbiology Data Center under the BioProject NMDC10018463. For each biological sample, the accession number as well as the corresponding link can be found in Supplementary Table 7. Sequences related to the parvoviruses identified in this study were previously submitted to GenBank under the submission numbers: OQ878961- OQ878963.

Additional information

Funding

This project was supported by the Chinese Ministry of Science and Technology (grant no. 2021YFC0863400), the Alliance of International Scientific Organizations (grant no. ANSO-CR-SP-2020-02), the Shanghai Municipal Science and Technology Major Project (grant no. 2019SHZDZX02), G4 funding from Institut Pasteur, Fondation Merieux and Chinese Academy of Sciences to G.W., and the International Affairs Department of the Institut Pasteur of Paris. A.T. is supported by the ANSO Scholarship for Young Talents. Y.K. is supported by the CAS-TWAS Fellowship for International Doctoral Students.