In April 2018, fecal samples of 31 birds and 41 rats were collected from various locations in the Cona County, Tibet Province (4,360 m above sea level; Supplementary Figure 1 and Supplementary Table 1). In July 2019, fecal samples of 25 birds and 50 small mammals were collected from various locations in Yushu and Nangqian Counties of Qinghai Province (3,890 m above sea level; Supplementary Figure 1 and Supplementary Table 1). Small mammals were captured in their natural habitat using mousetraps, and birds were coincidentally captured while catching pikas near pika holes (Zhu et al., 2021b). The wild animals were euthanized and dissected. Their fecal samples were collected and preserved in maintenance medium consisting of Hank’s balanced salt solution (pH 7.4) with penicillin G (100 U/ml) and streptomycin (50 μg/ml), and kept at −20°C while transfer to the laboratory, and stored at −80°C in the laboratory. The species of animals were identified using the mitochondrial cytochrome b (Cyt b) gene (Sorenson et al., 1999) or by morphological observation by experts. The sampling process was conducted by the local center for disease control and prevention (CDC) as part of the National Surveillance Program for Plague in Wildlife, and authorized by the Ethics Committee of National Institute for Communicable Disease Control and Prevention, China CDC (ICDC-2019012).

The process was mainly followed as previously reported (Wu et al., 2018). Briefly, each specimen was homogenized in phosphate buffer saline (PBS). Clear suspensions were obtained by centrifugation at 15,000 × g for 20 min and were filtered using a 0.22-μm polyvinylidene difluoride filter. The filtered supernatant was centrifuged at 300,000 × g for 2 h at 4°C. Pellets were re-suspended in PBS and digested using the RNase-Free DNase I Kit (Qiagen) at 37°C for 1 h. RNA was extracted using the QIAamp Viral RNA Mini Kit (Qiagen). RNA concentration and quality of each sample were determined using Qubit (Thermo Fisher) and 2100 Bioanalyzer (Agilent). The same animal classes from the same sampling county were divided into the same group (Supplementary Table 1), resulting in four groups. RNA of each sample of the same group was pooled in equal quantity.

To facilitate virus discovery, rRNA of each library was removed as previously descripted (Zhu et al., 2021a) using the Ribo-Zero Gold rRNA Removal Kit (Illumina). Libraries were constructed using the TruSeq Stranded Total RNA Library Prep Gold Kit (Illumina) according to manufacturer’s instructions. RNA was fragmented, and random hexamers were used to transcribe RNA into cDNA. The second strand of cDNA was obtained using the DNA polymerase I large fragment. The next steps included end repair, adapter ligation, purification, and fragment selection. The constructed libraries were sequenced using the Illumina HiSeq 2000 platform with 150 bp paired-end reads.

RNA viruses were detected according to an established metatranscriptomic pipeline (Shi et al., 2018). Raw reads were trimmed to remove the adapter and low-quality reads (<Q20) using Trimmomatic v0.32 (Bolger et al., 2014). The obtained high-quality reads were assembled de novo per pool using both Trinity v2.4.0 (Grabherr et al., 2011) and Megahit v1.1.2 (Li et al., 2015). The resulting contigs were first annotated using the database including all reference virus proteins downloaded from NCBI¹ by Diamond BLASTx with the e-values to 1e^–5 (Buchfink et al., 2015). In an attempt to identify highly divergent viral sequences, the assembled contigs were BLAST searched against the conserved domain database (CDD) v3.14 (Lu et al., 2020). Subsequently, the obtained viral contigs were verified and BLAST searched against both non-redundant protein and nucleotide databases, and contigs showing similarity to the host, plant, bacterial, and fungal sequences were eliminated. The resulting viral contigs were compared with their closely related members, and those with all viral proteins of the corresponding genus or family were retained for further analyses. The RNA-Seq by expectation–maximization algorithm (Li and Dewey, 2011) was used to quantify the abundance of contigs. Finally, viral contigs sharing <90% RNA-dependent RNA polymerase (RdRp) amino acid identity with any previously known virus were identified as new viruses (Shi et al., 2016).

Putative open reading frames (ORFs) of viral genomes were predicted using NCBI ORF Finder.² The putative function of viral protein was annotated using CDD.³ Genetic distance (p-dist) was estimated using MEGA X software (Kumar et al., 2018), and amino acid identity was calculated using BioAider v1.314 (Zhou et al., 2020). Codon usage preferences were estimated using the codon usage similarity index (COUSIN⁴).

Potential viral contigs and their closely related members were aligned by multiple alignment program using MAFFT v7 (Katoh and Standley, 2013). The best-fit substitution models were estimated using ModelFinder in IQ-TREE v2 (Minh et al., 2020). Phylogenetic trees were constructed using PhyML v3.0 based on the maximum-likelihood method (Guindon et al., 2010) with corresponding substitution models from ModelFinder and 1,000 bootstrap replicates. Trees were finally edited and visualized in interactive Tree of Life v1.0 (Letunic and Bork, 2016).

The probability of being able to infect humans was ranked using machine learning models with a 0.303 of the cutoff value following a recent report (Mollentze et al., 2021). Briefly, viral genomes were merged into a file in the FASTA sequence format, and the PredictNovel.R script was run.

Gaps between viral contigs with unassembled overlaps were filled by RT-PCR, which was performed using the PrimeScript™ One-Step RT-PCR Kit with specific primers (Supplementary Table 2) based on assembled sequences, and Sanger DNA sequencing. To confirm the assembly results, reads were mapped back to the viral sequences and aligned using Bowtie 2 (Langmead and Salzberg, 2012). To exclude the contigs belonging to expressed endogenous virus elements (EVEs), DNA was extracted from the corresponding samples and used for PCR amplification (Supplementary Table 2) using Taq DNA polymerase (TaKaRa). The sequence was eliminated if the PCR results were positive. In addition, these vial sequence (including host hits) were removed when performing BLASTn searches against the non-redundant nucleotide database.

To confirm the prevalence of vertebrate viruses, we designed specific primers (Supplementary Table 2) for viruses based on assembled genomes, and performed RT-PCR to screen corresponding viral sequence in individual samples. PCR products were subjected to gel purification and Sanger DNA sequencing. We did not submit the sequences from Sanger DNA sequencing to public databases.

A total of 184 complete or near-complete viral genomes that contain the complete RdRp domain were obtained (Figure 1 and Supplementary Table 3). Sequence comparisons indicated that 173 of them were divergent from previously known viruses, sharing <90% amino acid identity with known viruses. Coronaviruses detected in fecal samples have already been reported (Zhu et al., 2021c). Viruses of the family Picobirnaviridae accounted for >50% of the total number of viruses in the XZSCD and YSSCD libraries (Supplementary Figure 2B). All four libraries contained viruses from families Astroviridae, Iflaviridae, Partitiviridae, and Solemoviridae (Supplementary Figure 3).

To assess the amount of each viral RNA read, reads were mapped back to viral genomes. The viral family abundance among the four libraries showed marked differences (Supplementary Figure 2C). Viral families that revealed relatively high abundances were Nodaviridae, Chuviridae, Partitiviridae, Dicistroviridae and Solemoviridae in the XZNCD library; Solemoviridae, Chuviridae, Picornaviridae and Mitoviridae in the XZSCD library; Totiviridae in the YSNCD library; Nairoviridae in the YSSCD library (Supplementary Figure 2).

Despite the large number of viruses discovered, we mainly focused on the characteristics of vertebrate viruses, which signified birds or mammals as virus hosts. In addition to these insect, plant, and fungal virus families (Luo et al., 2018; Wu et al., 2018), viruses of families Rhabdoviridae, Phasmaviridae and Phenuiviridae appeared to be insect viruses, because they were grouped with invertebrate RNA viruses in the phylogenetic trees (Supplementary Figure 4). Thus, we mainly describe the characteristics and phylogenetic relationships of viruses within families Astroviridae, Hepeviridae, Nairoviridae, Picornaviridae and Picobirnaviridae.

Seven nearly complete astrovirus genomes (length ranging from 5,779 to 7,066 bp) were detected and assembled from the four libraries. Avastrovirus YSN01 and YSN02 were identified in one bird (Montifringilla taczanowskii) fecal sample from Qinghai Province. Mamastrovirus YSS01–YSS03 were detected in one (Apodemus peninsulae), two (A. peninsulae and Cricetulus kamensis), and two (A. peninsulae and Crocidura sp.) fecal samples from Qinghai Province, respectively. Mamastrovirus XZS01 and avastrovirus XZN01 were detected from six rat fecal samples (all positive samples were of Phaiomys leucurus) and one bird (Leucosticte brandti) fecal sample, respectively (Supplementary Table 4). Each of these detected astroviruses showed typical genome organization comprising a single-stranded positive RNA containing three ORFs. ORF1a and ORF1b encoded non-structural polyproteins (a protease and an RdRp), and ORF2 encoded the viral capsid precursor (Figure 2A). According to the phylogenetic analysis (Figure 2B) based on ORF2 protein sequences, avastrovirus YSN01, YSN02, and XZN01 were assigned to species Avastrovirus 5 of the genus Avastrovirus (Fernández-Correa et al., 2019), whereas, mamastrovirus YSS01 and YSS02 formed an independent clade. The genetic distances (p-dist) and amino acid identity of ORF2 between mamastrovirus YSS01, YSS02, and their closely related members were greater than the threshold value 0.741 and less than 75%, respectively, indicating that mamastrovirus YSS01 and YSS02 may represent a new species of genus Mamastrovirus. Mamastrovirus YSS03 and XZS01 were closely related to members of species Mamastrovirus 6 (Alves et al., 2018), but their p-dist (>0.741) and amino acid identity (<75%) values shared with closely related members indicate that they represent new species belonging to genus Mamastrovirus.

In the YSSCD library, one highly divergent virus (7,934 bp; MW826539) of the family Hepeviridae was identified that showed only 37.2% amino acid identity to known Hepeviridae sp. muf159hep1. Another nearly complete genome (6,940 bp; MW391926) was obtained and named Orthohepevirus C strain YS19. The virus was detected in two fecal samples of Ochotona curzoniae and A. peninsulae from Qinghai Province. Complete genome analysis showed that Orthohepevirus C strain YS19 was closely related to Orthohepevirus C isolates RdHEVAc86 (77.6% nucleotide identity) and RdHEVAc14 (77.3% nucleotide identity). The genome organization of YS19 shows four ORFs similar to those of RdHEVAc86 and RdHEVAc14 (Figure 3A and Supplementary Figure 5). The amino acid identity between YS19 and RdHEVAc86 and RdHEVAc14 ranged from 70 to 91% (Supplementary Table 5). A genome-based maximum-likelihood phylogenetic tree revealed that Orthohepevirus C strain YS19 formed an independent cluster along with RdHEVAc86 and RdHEVAc14 within the group of the species Orthohepevirus C (Figure 3B). Values of pairwise distances between YS19 and RdHEVAc86 and RdHEVAc14 were 0.068 and 0.069, respectively, which were below the threshold of 0.088 (Smith et al., 2016). These results indicate that YS19 belongs to the HEV-C3, which includes RdHEVAc86 and RdHEVAc14.

We identified an orthonairovirus in the family Nairoviridae from the YSSCD library, which was confirmed by PCR amplification using three paired primers (Supplementary Table 2). All three segments were found in only one shrew fecal sample from the Qinghai Province. The virus was named orthonairovirus YSS19 (YSV) under accession numbers MW391927–MW391929. The amino acid identity analysis indicated a close relationship of YSV with Erve orthonairovirus (ERVEV, 52.3–59.7%), Thiafora orthonairovirus (TFAV, 51.2–59.4%), and Crimean-Congo hemorrhagic fever orthonairovirus (CCHFV, 47.6–48.8%). Phylogenetic analysis based on the putative protein sequence of the L segment showed that YSV grouped with ERVEV and TFAV, which were closely related to CCHFV (Figure 3C). In addition, phylogenetic analysis based on putative protein sequences of M and S segments also showed similar results (Supplementary Figures 6, 7).

Codon usage preferences of YSV, ERVEV, TFAV, and CCHFV showed a highly similar pattern (Supplementary Figure 8). The glycoprotein precursor of YSV includes one unique proteolytic site and another identical to the proteolytic site of CCHFV and Nairobi sheep disease virus (NSDV) groups. The ovarian tumor domain (Frias-Staheli et al., 2007; Dilcher et al., 2012) of YSV (116 amino acids) shared 60.2%, 55.9%, and 45.8% amino acid identity with domains of ERVEV, TFAV, and CCHFV, respectively (Supplementary Figure 9).

Eight nearly complete picorna-like virus genomes were obtained (Figure 4A). Three viruses [Picornaviridae sp. XZS01 and XZS02, and rat sapelovirus (SPV) XZS03] were detected in rat fecal samples from Tibet Province, three viruses [rat mosavirus YSS01, rat SPV YSS02, and hepatovirus A (HAV) isolate YSS03] were identified in rat fecal samples from Qinghai Province, and two viruses (bird SPV YSN01 and Picornaviridae sp. YSN02) were found in bird fecal samples from Qinghai Province. A virus with <33% and 36% amino acid identity in P1 protein and 2C + 3 CD proteins may belong to a new genus of the family Picornaviridae (Zell et al., 2017). Rat SPV YSS02 and XZS03 shared 55.0% and 61.4% polyprotein amino acid identity with known members of genus Sapelovirus, indicating that they represent putative new species of that genus (Figure 4B). Polyprotein of rat mosavirus YSS01 shared 61.4% amino acid identity to that of Mosavirus A2. Picornaviridae sp. YSN02 shared 50.9% polyprotein amino acid identity with known Picornaviridae sp. sru051pic1. Bird SPV YSN01, Picornaviridae sp. XZS01, and Picornaviridae sp. XZS02 shared <40.5%, <38.7%, and <31.4% polyprotein amino acid identity with all known viruses of the family Picornaviridae, respectively, indicating that these three viruses may belong to putative new species or genera in family Picornaviridae.

HAV YSS03 (7,575 bp; MW391925) was identified in two samples from A. peninsulae; its genome included a 719-nucleotide 5′-untranslated region, 6756-nucleotide ORF-encoding polyprotein precursor (2,251 amino acids), and 100-nucleotide 3′-untranslated region. Sequence comparison of the HAV YSS03 genome showed 98.1% and 97.8% nucleotide similarities with Marmota himalayana hepatovirus (MHHAV) 2ID and 3ID, respectively. Interestingly, HAV YSS03 also shared 99.6% and 99.3% amino acid identity with MHHAV 2ID and 3ID, respectively. The potential polyprotein cleavage sites were identical between HAV YSS03 and MHHAV (Yu et al., 2016). Phylogenetic analysis based on genomic sequences (Supplementary Figure 10) revealed that HAV YSS03 clustered with MHHAV forming an independent branch within genus Hepatovirus. These results indicated that a nearly identical hepatovirus was present in sympatric yet genetically different hosts.

The family Picobirnaviridae includes only one genus, Picobirnavirus, which is usually found in vertebrate fecal samples (Luo et al., 2018). Viruses possess bi-segmented double-stranded RNA genomes (i.e., segments 1 and 2) and are rarely unsegmented (Luo et al., 2018; Delmas et al., 2019). A total of 29 segment 1 (containing the RdRp domain) and 39 segment 2 (capsid) sequences were identified in rat fecal samples from Tibet Province with the addition of two unsegmented picobirnaviruses (PBVs). Moreover, a total of 28 segment 1 and nine segment 2 sequences were identified in rat fecal samples from Qinghai Province. The amino acid identity of RdRp with that of known PBVs was <83.6%, including only two RdRp sequences showing >80% identity. To investigate the evolutionary position of these PBVs, phylogenetic trees based on RdRp and capsid amino acid sequences were constructed (Figure 5A), and they revealed that newly discovered capsid sequences belong to eight clusters (C1–C8). In addition to known genotypes (GI–GVI), these newly discovered RdRp sequences also formed two new genogroups (GVII and GVIII). The six members of the GVII genogroup (Figure 5A) were detected in fecal samples from Tibet Province, sharing 23.6–43.1% amino acid identity with known PBVs. The GVIII genogroup contained only one virus (Rat PBV XZ01) with 48.1% amino acid identity with its closest relative. Interestingly, two unsegmented PBVs, named rat PBV XZ03 (3,494 bp) and rat PBV XZ04 (4,494 bp), were identified and had complete ORFs encoding the RdRp and capsid regions. The assortment types (Figure 5B) of rat PBV XZ03 (C2: GI) and rat PBV XZ04 (C8: GV) were different from those previously reported (Luo et al., 2018). Rat PBV XZ03 was detected in two rat (Phaiomys leucurus) fecal samples, and rat PBV XZ04 was detected in one fecal sample.

Zoonotic risk (i.e., the probability of being able to infect humans) of newly discovered vertebrate viruses was evaluated as previously reported (Mollentze et al., 2021). Ten viruses (avastrovirus XZN01, avastrovirus YSN02, mamastrovirus YSS02, mamastrovirus YSS01, rat SPV YSS02, bird SPV YSN01, Picornaviridae sp. XZS02, HAV YSS03, YSV, and Orthohepevirus C strain YS19) were high priority (Figure 6). Mamastrovirus XZS01, rat SPV XZS03, and Tibet bird virus 2 were ranked as low priority (Figure 6).

In addition to the viruses described above, a large number of invertebrate RNA viruses were detected. These viruses belonged to the families Partitiviridae (n = 23), Solemoviridae (n = 21), Iflaviridae (n = 8), Dicistroviridae (n = 7), Nodaviridae (n = 7) and Totiviridae (n = 6) (Figure 1). Also, 12 of these viruses were distantly related to previously known viruses and each other. The RdRp domain of these viruses shared 26.3–37.6% amino acid identity with corresponding sequences of their closest relatives (Table 1), indicating that these 12 viruses may belong to potentially new genera or families.