Molecular detection and characterization of three novel parvoviruses belonging to two different subfamilies in zoo birds

Birds carry a large number of viruses that may cause diseases in animals or human. At present, virome of zoo birds are limited. In this study, using viral metagenomics method, we investigated the feces virome of zoo birds collected from a zoo of Nanjing, Jiangsu Province, China. Three novel parvoviruses were obtained and characterized. The genome of the three viruses are 5,909 bp, 4,411 bp and 4,233 bp in length respectively which encoded four or five ORFs. Phylogenetic analysis indicated that these three novel parvoviruses clustered with other strains formed three different clades. Pairwise comparison of NS1 amino acid sequences showed that Bir-01–1 shared 44.30%−74.92% aa sequence identity with other parvoviruses belonging to the genus Aveparvovirus, while Bir-03-1 and Bir-04–1 had lower than 66.87% and 53.09% aa sequence identity with other parvoviruses belonging to the genus Chaphamaparvovirus. These three viruses were identified as three novel species of the genus Aveparvovirus and Chaphamaparvovirus respectively basing on the species demarcation criteria of parvovirus. Our findings broaden the knowledge of the genetic diversity of parvovirus and provide epidemiological data for the outbreak of potential bird’s parvovirus disease.

known if they are associated with avian diseases [4]. The genus Aveparvovirus currently includes four species -Aveparvovirus columbid1, Aveparvovirus galliform1, Aveparvovirus gruiform1, and Aveparvovirus passeriform1 -but many other aveparvovirus-like viruses have been detected recently in various wild birds [5,6], suggesting that new species may need to be added to the genus Aveparvovirus.
In 2018, cloacal swabs of 58 birds of 14 species from three orders (Gruiformes, Galliformes, and Psittaciformes) were collected from a zoo in the city of Nanjing in Jiangsu province (Supplementary Table S1). Cloacal swabs were collected by inserting disposable absorbent cotton swabs approximately 3 cm into the cloaca, and these were stored on dry ice prior to shipping to the laboratory. Sample preparation for metagenomic sequencing was performed as described previously. The swab tips were immersed into 500 µL of Dulbecco's phosphate-buffered saline (DPBS) and vortexed vigorously for 10 minutes. The supernatants were then collected in new 1.5-ml Eppendorf tubes, centrifuged at 15,000 × g for 10 minutes, and stored at -80°C. For viral metagenomic sequencing, four libraries were constructed. A 50-µL aliquot of the supernatant was removed from each tube and combined to form four pools (700-750 µL of supernatant per pool). Sample pools were mixed and centrifuged at 12,000 × g for 20 min at 4°C. The supernatant was passed through a 0.45-µm filter (Merck Millipore, MA, USA) to remove eukaryotic and bacterial cell-sized particles. Filtrates were then digested with DNase and RNase at 37°C for 60 min. Total nucleic acid was then extracted using a QIAamp MinElute Virus Spin Kit (QIAGEN, HQ, Germany) according to the manufacturer's protocol.
The enriched viral nucleic acids from each pool were used as templates for reverse transcription reactions, using reverse transcriptase and 6-base random primers, followed by a single round of DNA synthesis using Klenow fragment polymerase (New England BioLabs, MA, USA). Four libraries were constructed using a Nextera XT DNA Sample Preparation Kit (Illumina, CA, USA) and sequenced on an Illumina Miseq platform with 250 base-paired ends with dual barcoding for each pool. For bioinformatic analysis, barcodes of raw sequences were removed using the vendor software from Illumina. 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 to 55 were identical, and only one random copy was kept. Lowsequencing-quality tails were trimmed, using Phred quality score 10 as the threshold. The total number of reads in each library is shown in Table 1. The adapter sequence was removed using NCBI VecScreen (http://www.ncbi.nlm.nih. gov/tools/vecscreen/) with special parameters for adapter removal. Bacterial reads were subtracted by mapping to the bacterial BLAST NT database using bowtie2 v2.2.4. The clean reads were assembled de novo using SOAPdenovo2 r240 with a kmer size of 63 and default settings. The assembled contigs, along with singlets, were aligned with an inhouse viral proteome database using BLASTx (v.2.2.7) with an E-value cutoff of < 10 − 5 . The viral BLASTx database was compiled using a fasta file of the NCBI virus reference proteome (https://ftp.ncbi.nih.gov/refseq/release/viral/).
A total of 9,099,916 reads were obtained from four libraries. Contigs and singlets were matched against a customized viral proteome database using BLASTx with an E value cutoff of < 10 − 5 . Of those, 6,829 reads were assigned to the family Parvoviridae. Three complete genome sequences of parvoviruses were obtained by a combination of assembly and PCR to bridge sequence gaps, and these were named Bir-01-1 (mean coverage, 1,134), Bir-03-1 (mean coverage, 25), and Bir-04-1 (mean coverage, 126.7). The three parvovirus sequences were deposited in the GenBank database with accession numbers OP894050 to OP894052. The raw sequence reads from the metagenomic library were deposited in the Short Read Archive of the GenBank database (Supplementary Table S1).
The genomes of the three viruses are 5,909, 4,411, and 4,233 nt in length, respectively. ORFs were predicted using a combination of Geneious 11.1.2 software and NCBI ORFfinder. The results showed that Bir-01-1 and Bir-04-1 both contain five open reading frames (ORFs), while Bir-03-1 has four ORFs (Fig. 1A). Bir-01-1 encodes one nonstructural protein (NS1) and two structural proteins (VP1 and VP2), while Bir-03-1 and Bir-04-1 encode three nonstructural proteins (NS1, NS2, and NS3) and one structural protein (VP1). The NS1 protein of these three novel parvoviruses are 640 aa, 633 aa, and 654 aa in length, respectively, and possess two conserved replication initiator motifs (HuH and Yxx(x)K) and a conserved NTPase/ helicase motif (Walker A, B, B', and C) (Fig. 1B) as in other parvoviruses [1]. BLASTp analysis based on the NS1 protein showed that Bir-01-1 had the highest aa sequence identity of 99.22% to isolate bpk075par01 (GenBank no. MT138216) from China but had only 44.30%-74.92% aa sequence identity to other parvoviruses belonging to the  Fig. S1). Based on the species demarcation criteria proposed by Pénzes and co-workers, viruses can be considered members of the same species if the NS1 proteins share more than 85% aa sequence identity, while NS1 proteins of members of the same genus should share at least 35-40% aa sequence identity with a coverage of > 80% between any two members [1]. According to these criteria, Bir-01-1 belongs to a novel species of the genus Aveparvovirus, while Bir-03-1 and Bir-04-1 belong to two different species of the genus Chaphamaparvovirus. The VP1 protein of these three strains are 675 aa, 529 aa, and 532 aa in length, respectively, and lack the conserved phospholipase A2 (PLA2) motif that is associated with virus infectivity, suggesting that they might use a different mechanism for virus entry and release [7]. BLASTp analysis based on the VP1 protein showed that Bir-01-1 shared the highest aa sequence identity of 87.41% with the isolate bpk075par01 (GenBank no. MT138216) from China. A glycine-rich sequence ( 163 GAGGGGGGGGVG 174 ) resembling those of other members of the genus Aveparvovirus [8] was found in the VP1 protein of Bir-01-1. VP1 of Bir-03-1 and Bir-04-1 had the highest aa sequence identity of 76.50% and 66.34% to isolate par083par07 (GenBank no. MW046512) from China and isolate yc-9 (GenBank no. KY312548) from a fecal sample from a red-crowned crane in China, respectively. No glycine-rich sequences were found in the VP1 proteins of Bir-03-1 and Bir-04-1. Phylogenetic analysis was performed based on the deduced amino acid sequences of the NS and VP proteins of the parvovirus with the closest BLASTx matches in Gen-Bank and representative members of related viral species or genera. Phylogenetic trees with 1000 bootstrap resamplings of the alignment data sets were generated using the neighbor-joining method in MEGA 11. Phylogenetic trees were also constructed using MrBayes v3.2.7. The Markov chain was run for a maximum of one million generations, in which every 50 generations were sampled and the first 25% of Markov chain Monte Carlo samples were discarded as burn-in (Supplementary Fig. S2). The results showed that samples tested (10/58) were positive for Bir-01-1, which was detected in samples from blue peacock, red-waisted golden pheasant, Reeve's pheasant, silver pheasant, cockatoo, grey crane, and white-naped crane; 8.62% of fecal samples (5/58) were positive for Bir-03-1, which was detected in samples from blue-and-yellow macaw, gray parrot, and red-crowned crane; and 17.24% of fecal samples (10/58) were positive for Bir-04-1, which was detected in samples from white crane, gray parrot, and red-crowned crane. It has been reported that aveparvoviruses were widely distributed in wild birds and poultry, with a prevalence rate of 12.1-78% [9]. In this study, aveparvoviruses were detected in birds with a positive rate of 17.2%, which is consistent with previous studies. It has been reported that chaphamaparvoviruses may be host-specific. The results of this study showed that four of the samples that were positive for Bir-03-1 were from members of the order Psittaciformes, but one was from a member of the order Gruiformes. Similarly, 11 samples that were positive for Bir-04-1 were from birds of the order Gruiformes. Bir04-1 and the parvovirus found in Chinese red-crowned cranes (Gruiformes) clustered in the same branch [10]. This provides further evidence of the host specificity of chaphamaparvoviruses. Chaphamaparvoviruses may be pathogenic to mammals, since they have been found in inclusion bodies of kidney cells of mice with NS1 of Bir-01-1 is phylogenetically related to pigeon parvovirus A (KC876004), and it clustered with strains detected in wild birds (pigeon, parrot, red-crowned crane) or poultry (chicken and turkey) in China and the USA to form a clade within the genus Aveparvovirus, subfamily Parvovirinae. NS1 of Bir-03-1 and Bir-04-1 clustered with strains in the newly established subfamily Hamaparvovirinae. NS1 of Bir-03-1 was found to be most closely related to yc-9 (KY312548), which was detected in a red-crowned crane in China, and NS1 of Bir-04-1 clustered with BR_DF10 (MN175612), which was detected in a neognath bird in Brazil ( Fig. 2A). The phylogenetic tree based on VP1 sequences was consistent with the NS1-based tree (Fig. 2B).
To investigate the prevalence of parvoviruses in zoo birds, sets of primers were designed to amplify the NS1 region of these three parvoviruses ( Table 1). For screening of Bir-01-1, the forward primer 5-CCCCTAACAGGCA-CAGAAGG-3 and the reverse primer 5-GGTTGATTT-GTTGCGGGACC3 were used, for screening of Bir-03-1, the forward primer 5-ATCTGATGCTTGCAGACGCT-3 and the reverse primer 5-ATTCACCACCTGTGCATCGT-3 were used, and for screening of Bir-04-1, the forward primer 5-CGAATGCCGTTATTGCCGAG-3 and the reverse primer 5-TCCCCAACAGTCCCATTGTG-3 were used ( Table 1). The results showed that 17.24% of the fecal kidney disease, and they have also detected in the feces of dogs with diarrhea [11,12]. However, further animal testing or histopathological examination will be needed to determine whether they can cause disease in birds. Coinfection with Bir-01-1 and Bir-04-1 was detected in one sample from grey crane. Because no obvious clinical symptoms were observed in the zoo birds during sampling, the possible relationship between these three parvoviruses and disease needs to be clarified. In this study, in addition to the three parvoviruses, we detected other viruses, such as duck hepatitis A viruses and CRESS DNA viruses. Here however, we specifically focused on the parvoviruses due to the high proportion of parvovirus sequence reads in the four libraries and the potential pathogenicity of these viruses to birds. Further analysis and research on the other viruses will be conducted in subsequent studies.
In summary, we identified three novel parvoviruses in different zoo birds and determined their complete genome sequences. The NS1 proteins of these three novel parvoviruses possess conserved replication initiator and NTPase/ helicase motifs. Phylogenetic analysis based on NS1 and VP1 sequences indicated that Bir-01-1 clustered with other strains belonging to the genus Aveparvovirus to form a clade, while Bir-03-1 and Bir-04-1 clustered with members of the newly established subfamily Hamaparvovirinae to form two different clades. Based on the latest demarcation criteria for parvoviruses, Bir-01-1 belongs to a novel species in the genus Aveparvovirus, while Bir-03-1 and Bir-04-1 belong to two different species in the genus Chaphamaparvovirus. The epidemiological results suggest that these three viruses are prevalent in zoo birds.