From May to July 2016, fresh fecal specimens of seven healthy adult wild giant pandas were collected using disposable absorbent cotton swabs in the southern slope area of the middle section of the Qinling Mountains at an altitude of about 2,500–3,000 m of Mainland China. Fecal samples were collected using disposable materials. All specimens were stored in sterilized covered containers and transported over dry ice. Prior to viral metagenomic analysis, 10 g of each fecal sample was suspended using 5 ml Dulbecco's phosphate-cushioned saline (DPBS) and vortexed for 5 min vigorously, then incubated at 4°C for 30 min. After centrifugation (10 min, 15,000 g), 1 ml of each supernatant was collected in 1.5-ml centrifuge tubes and stored at −80°C for later use (15). Since these giant pandas are wild, we do not know their exact age. Their approximate age, physical condition, and freshness of the feces were inferred by local experienced forest workers based on the characteristics of feces.

A total of 500 μl of supernatant was pipetted from the seven samples (71.43 μl per sample was a pipetted from) and collected into a new 1.5-ml tube. The samples were centrifuged at 4°C for 5 min at 12,000 g and then filtered through a 0.45-μm filter (Merck Millipore, MA, USA) to remove bacterial cell-sized and eukaryotic particles. Filtrates were then treated with DNase (Turbo DNase from Ambion, Thermo Fisher, Waltham, MA, USA; Baseline-ZERO from Epicentre, Charlotte, NC, USA; and Benzonase from Novagen, Darmstadt, Germany) and RNase (Promega, Madison, WI, USA) to digest unprotected nucleic acids at 37°C for 60 min (16). Total nucleic acids were then extracted using QIAamp MinElute Virus Spin Kit (Qiagen, HQ, Germany) according to the manufacturer's protocol. These nucleic acid samples containing DNA and RNA viral sequences were subjected to reverse transcription reactions with SuperScript III reverse transcriptase (Invitrogen, Carlsbad, CA, USA) and 100 pmol of a random hexamer primer, followed by a single round of DNA synthesis using Klenow fragment polymerase (New England Biolabs, Ipswich, MA, USA). One library was constructed using the Nextera XT DNA Sample Preparation Kit (Illumina, CA, USA) and sequenced on the MiSeq Illumina platform with 250 base paired ends with dual barcoding. Paired end reads of 250 bp generated by MiSeq were debarcoded using vendor software from Illumina for bioinformatics analysis. An in-house analysis pipeline running on a 32-node Linux cluster was used to process the data. Reads were considered duplicates if bases 5–55 were identical and only one random copy of duplicates was kept. Clonal reads were removed, and low sequencing quality tails were trimmed using Phred quality score 20 as the threshold. Adaptors were trimmed according to the default parameters of VecScreen, which is the National Center for Biotechnology Information (NCBI) BLASTn with specialized parameters designed for adapter removal. The cleaned reads were de novo assembled within each barcode using the ENSEMBLE assembler (17). All generated contigs and singlet reads were then matched against a customized viral proteome database using BLASTx with an E-value cutoff of <10⁻⁵, where the virus BLASTx database was compiled using the NCBI virus reference proteome (https://ftp.ncbi.nih.gov/refseq/release/viral/) to which viral protein sequences from NCBI nr fasta file (based on annotation taxonomy in Virus Kingdom) were added. Contigs without significant BLASTx identity to the viral proteome database are searched against viral protein families in the vFam database (18) using HMMER3 to detect remote viral protein similarities (19–21).

To assemble the potential viral genome of interest, the contigs which showed significant BLASTx identity to related viruses were selected in Geneious software version 11.1.2 (22). The contigs with consensus sequence length >500 bp were subjected to further analysis where the individual contig was used as reference for mapping to the raw data of its original barcode using the Low Sensitivity/Fastest parameter in Geneious 11.1.2. Those prolonged contigs with sequence length >1,500 bp and the major non-structural protein and the structural protein, as well as some contigs which only had a non-structural protein, were subjected to full characterization. Putative viral open reading frames (ORFs) were predicted by Geneious v11.1.2 with built-in parameters (minimum size: 300; genetic code: standard; start codons: ATG) (22) and were further checked through comparing to related viruses. The annotations of these ORFs were based on comparisons to the Conserved Domain Database. Putative exon and intron were predicted by Netgenes2 at http://www.cbs.dtu.dk/services/NetGene2/.

PCR screening was performed for the novel viruses of Picornavirales in the fecal samples to facilitate the investigation of prevalence. Sequences and characteristics of the primers used in this study are shown in Supplementary Table S2.

To infer phylogenetic relationships, nucleotide sequences of reference strains belonging to different groups of corresponding viruses were downloaded from the NCBI GenBank database. Related nucleotide sequences were aligned using the alignment program implemented in the CLC Genomics Workbench 10.0, and the resulting alignment was further optimized using MUSCLE in MEGA v7.0 (23) and MAFFT v7.3.1 employing the E-INS-I algorithm (24). Sites containing more than 50% gaps were temporarily removed from alignments. Bayesian inference trees were then constructed using MrBayes v3.2 (25). The Markov chain was run for a maximum of 1 million generations, in which every 50 generations were sampled and the first 25% of Markov chain Monte Carlo (mcmc) samples were discarded as burn-in. Maximum likelihood trees were also constructed to confirm all the Bayesian inference trees using software MEGA v7.0 (23).

Standard precautions were used for all procedures to prevent cross-sample contamination and nucleic acid degradation. Aerosol filter pipette tips were used to reduce the possibility of sample cross-contamination, and all the materials (including microcentrifuge tubes, pipette tips, etc.) which directly contacted with nucleic acid samples were RNase and DNase free. Nucleic acid samples were dissolved in DEPC-treated water, and RNase inhibitors were added. As a blank control, sterile ddH2O (Sangon, Shanghai, China) was simultaneously prepared and further processed in the same condition. During quality inspection using agarose gel electrophoresis and Agilent Bioanalyzer 2100, no detectable DNA existed in the control library. During sequencing of the Illumina MiSeq platform, the control library generated a small number of sequence reads. BLASTx searching based on the total reads in the control library revealed no viral sequences.

The raw sequence reads from the metagenomic library were deposited in the Short Read Archive (SRA) of the GenBank database with accession no. SRR11348859.

All sequences obtained in this study have been deposited in the GenBank database under accession numbers MZ956587-MZ956589, MZ956591, and MZ956593-MZ956598.

The library generated a total of 393,698 raw sequence reads on the Illumina MiSeq platform; the number of unique reads after removing duplicates was 131,245. A total of 1,411 contigs were generated, and the average length of contigs was 710 bp. After bioinformatics analysis and mapping for Megan6 software visualization, a total of 35,800 sequence reads had the best matches with viral proteins (Supplementary Table S1), accounting for 9.09% of the total raw sequence reads. About 40 viral families were detected; the most abundant virus family was Drexlerviridae (51.29% of the total analyzed virus reads), followed by Siphoviridae (17.93%), Ackermannviridae (5.79%), Genomoviridae (5.35%), Microviridae (5.07%), Iflaviridae (4.31%), Dicistroviridae (3.88%), and Myoviridae (2.20%) (Figure 1 and Supplementary Table S1). In addition, there are 26,507 reads that have not been clearly assigned (E value >10⁻⁵), and potentially novel viruses are likely to exist in these reads.

In this study, six different RNA virus sequences (temporarily named PPLV1-PPLV6) related to Picornavirales were obtained using the assemble sequence program in Geneious 11.1.2; genome sizes of these six PPLVs were 9,653, 9,625, 9,670, 8,747, 4,694, and 8,638 bp, with GC contents of 41.3, 36.0, 34.3, 41.5, 34.9, and 38.3%, respectively, and their amino acid sequence identities with the best matches in GenBank ranged from 27.19 to 96.68%. The genomes of PPLV1, PPLV3, and PPLV5 all have a large ORF. PPLV1 and PPLV3 contain the Capsid domain, RNA helicase domain, and RdRp region, while PPLV5 has only one RdRp region. The genomes of PPLV2, PPLV4, and PPLV6 contain two ORFs, ORF1 and ORF2, respectively. In addition to the common RNA helicase domain, ORF1 of these three sequences also contains other different structures; for example, PPLV2 has a GTPase region, and PPLV4 has a 3C protease region, but their ORF2 all contain Capsid domain (Figure 2A). Primers were designed based on the RdRp sequence. PCR screening results showed that one sample was positive for PPLV1 and PPLV3 at the same time, and the other four viruses existed in the other four samples, respectively (Supplementary Table S2).

Phylogenetic analysis based on RdRp amino acid sequences indicated that PPLV1 clustered in the clades between the genus Iflavirus and Picorna-like virus, which could be considered a putative novel genus obviously. Similarly, PPLV3, as a separate branch, clustered with the genus Iflavirus. In addition, PPLV4, as an independent clade, was close to the Picorna-like virus. PPLV6 was clustered with the previously reported unclassified Dicistroviridae identified in the bird's intestines and has a 96.68% amino acid sequence identity. PPLV2 and PPLV5 shared 30.92 and 27.19% amino acid sequence identity with the best matches in the GenBank database, respectively. Furthermore, PPLV5 was classified outside the currently well-defined families, which makes it possible for the existence of a novel virus family (Figure 2B).

In this study, we have identified four complete genomovirus-related sequences through assembly and mapping, namely, Pgenomo1-4. The whole genome sizes of Pgenomo1-4 were 2,188, 2,213, 2,198, and 2,142 bp, respectively (Figure 3A). The stem-loop structure was found between the 5′ ends of the Rep protein and Cap protein of the genomes of all the four viruses (Figure 3B), all of which contained a conserved non-amer of “TARTRTTK” (Figure 3C). BLASTp searching in the GenBank database with the predicted Rep protein sequence indicated that the best match of Pgenomo1 was Luscinia calliope Genomoviridae sp. (QTE03633) from wild bird feces, sharing 90.69% amino acid sequence identity. The Rep protein of Pgenomo2 was closely related to Turdus pallidus Genomoviridae sp. (QTE03681), sharing 67.92% amino acid identity with it, but the genome-wide nucleotide identity was 81.70%; based on the genome-wide pairwise sequence identity threshold of 78%, this virus cannot be considered as a new species (26). The amino acid sequence identity between Pgenomo3 and its best matching Tarsiger cyanurus Genomoviridae sp. (QTZ83255) was 81.56%. Similarly, the best matches of Pgenomo4 was Phylloscopus Genomoviridae sp. (QVW56506) from wild birds, sharing 91.32% amino acid sequence identity (Figure 3D).