In 2020, with the assistance of animal experts, 10 healthy bar-headed geese were successfully captured using a cannon net in Saya County, Shigatse City, Tibet, China. Subsequently, fecal samples were collected using disposable, absorbent cotton swabs and shipped to our laboratory on dry ice. These samples were resuspended in 2 milliliters of Dulbecco’s Phosphate-Buffered Saline (DPBS) and vigorously vortexed for 10 min, followed by three freeze–thaw cycles. After centrifugation (15,000 × g, 4°C, 10 min), the supernatant of each sample was collected in a new 1.5 milliliter centrifuge tube and stored at-80°C for further use. Ethical approvals were given by the Ethics Committee of Jiangsu University with the reference number 2018ujs18023. Sample collection was performed in accordance with the Wildlife Protection Law of the People’s Republic of China. After sampling, the bar-headed geese were released without any damage.

A total of 500 μL of fecal suspension (50 μL of fecal supernatant from each fecal sample) was pooled together and filtered through a 0.45 μm filter (Merck Millipore, Billerica, MA, USA) to remove bacterial and eukaryotic cell-sized particles. The filtrates were then treated with a mixture of nuclease enzymes to digest the unprotected nucleic acids at 37°C for 90 min. Viral RNA and DNA were extracted by using the QIAamp MinElute Virus Spin Kit (Qiagen, Hilden, NRW, Germany) according to the manufacturer’s improved protocol. The concentrations of DNA and RNA were calculated using the Qubit 4 (Invitrogen, Carlsbad, CA, USA) nucleic acid concentration sequencer. The RNA and DNA were stored at-80°C for further use.

The viral nucleic acid pool containing DNA and RNA viral sequences was subjected to RT reactions with SuperScript III reverse transcriptase (Invitrogen, CA, USA) using 100 pmol of a random hexamer primer. The RT reaction conditions are 25°C for 10 min, 50°C for 60 min, 85°C for 5 min, and 95°C for 2 min. Then, the reaction products were quickly removed and placed on ice for >2 min. The klenow enzyme (New England Biolabs, MA, USA) was used to generate the complementary chain of cDNA. The Klenow reaction conditions are 37°C for 60 min, 75°C for 20 min. Next, libraries were constructed using the Nextera XT DNA Sample Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s protocol. A brief summary of the process includes adding connector primers and conducting 15 cycles of limited amplification. The library sequencing was completed by the Personalbio company using the NovaSeq Illumina platform with 250 base-paired ends with dual barcoding for each pool (12).

For bioinformatics analysis, the paired-end reads of 250 bp generated by NovaSeq were debarcoded using the vendor software from Illumina, which was processed using an in-house analysis pipeline running on a 32-node Linux cluster. Low-sequencing-quality tails were trimmed using a Phred quality score threshold of 10. Adaptors were trimmed using the default parameters of VecScreen, which is an NCBI BLASTn program with specialized parameters designed for adaptor removal. Bacterial reads were subtracted by mapping them to bacterial nucleotide sequences from the BLAST NT database using Bowtie2 v2.2.4. The cleaned reads were then de novo assembled by SOAPdenovo2 version r240 using a Kmer size of 63 with default settings. The assembled contigs, along with singlets, were aligned to an in-house viral proteome database using BLASTx (v.2.2.7) with an E-value cutoff of <10⁻⁵. The candidate viral hits were compared to an in-house nonvirus nonredundant (NVNR) protein database to remove false-positive viral hits. The NVNR database was compiled using nonviral protein sequences extracted from the NCBI NR Fasta file and was based on a taxonomy annotation excluding the virus kingdom.

The analysis of evolutionary relationships was carried out using amino acid sequences predicted from the genomic data, with reference to the closest viral relatives determined by the best BLASTx hit and representative members of related viral species or genera. Sequence alignment was conducted with Clustal W in MEGA-X using the default settings (13). Phylogenetic trees were constructed using MrBayes v3.2.7. The parameter “prset aamodelpr = mixed” was employed to enable the program to use the ten built-in amino acid models. The maximum number of generations was set to be ten million, and sampling occurred at every 50 generations, with the first 25% of Markov chain Monte Carlo (MCMC) samples being discarded during burn-in. Convergence was confirmed when the standard deviation of split frequencies was below 0.01. Bootstrap values were assigned to each node.

The pairwise comparison of viral amino acid sequences was conducted using SDTv1.2 software with default settings. Putative open reading frames (ORFs) in the genome were predicted using Geneious 11.1.2 software and the NCBI ORF finder.

The complete viral genome sequences identified in this study were deposited in GenBank under the accession numbers PQ220257-PQ220293. The raw sequence reads from the metagenomic library were deposited in the Sequence Read Archive of the GenBank database under accession number SRR29463689.

This study generated a total of 2,642,510 raw sequences using the Illumina NovaSeq platform. After bioinformatics analysis, 107,069 sequence reads were identified as viral sequences, accounting for 4.05% of the total raw data reads. Further analysis determined the proportion of these viral reads among different virus families. The result showed that the Picornaviridae family accounted for the highest proportion at 62%, followed by Parvoviridae at 19%, Circoviridae at 17%, and both Caliciviridae and Astroviridae at 1%, respectively (Figure 1). Because picornaviruses were analyzed in our previous study, this study focuses on other newly discovered viruses.

Astroviruses are a large group of non-enveloped, single-stranded positive-sense RNA viruses with spherical virions approximately 28–30 nanometers in diameter. First identified in 1975 in the feces of children with diarrhea, astroviruses have since been found to infect dozens of animal species. The Astroviridae family is divided into two genera: Mamastrovirus, which infects mammals, and Avastrovirus, which infects birds.

In this study, 519 viral sequence reads matching parvoviruses were detected from the bar-headed goose sample library. Using Geneious 11.1.2 for assembly, four partial astrovirus genomes were obtained and designated as ASTV1, ASTV4, ASTV5, and ASTV19. The lengths of these genomes are 5,679 nt, 2,370 nt, 2,162 nt, and 2,603 nt, respectively. ASTV1 encodes two complete ORFs (ORF1a and ORF1b) and a truncated capsid protein at the 3′ end. ASTV4, ASTV5, and ASTV19 have different lengths of coding regions: ASTV4 encodes incomplete ORF1a and capsid protein, while ASTV5 and ASTV19 encode parts of ORF1b and the capsid protein region (Figure 2A). BLASTx analysis of the sequences in NCBI revealed varying levels of sequence identity. ASTV1 showed the highest sequence identity (63.45%) with a strain isolated from mute swan feces (MW588064). ASTV4 and ASTV5 had 70.67 and 70.62% sequence identity, respectively, with a strain from Pacific black crow feces in Australia (MT894397). ASTV19 exhibited the highest identity (99.67%) with a strain previously isolated in our laboratory from wild bird feces (MT138004).

A phylogenetic tree was constructed based on the shared ORF1b of ASTV1, ASTV4, ASTV5, and ASTV19, along with representative strains from different genera of the Astroviridae family. As shown in Figure 2B, these four viruses cluster into two distinct branches. ASTV4 clusters with the Australian duck virus detected in 2020 (MT894397). Meanwhile, ASTV1, ASTV5, and ASTV19 cluster together with two viruses isolated from turkeys in the United States in 2012 (JQ307841 and JQ307840), as well as a virus detected in birds in China (MT137998). Additionally, a conserved TEEEY protein domain was identified in the ORF1a of ASTV1 and ASTV4 like other avastroviruses (Figure 2C).

The Caliciviridae family consists of a large group of non-enveloped, single-stranded positive-sense RNA viruses. Their genomes range from 6.4 to 8.5 kb in length and contain multiple ORFs encoding both non-structural and structural proteins. The virus particles exhibit icosahedral symmetry and have diameters ranging from 27 to 40 nm. This family is divided into 11 genera, including Lagovirus, Norovirus, Nebovirus, Recovirus, Sapovirus, Valovirus, Vesivirus, Bavovirus, Nacovirus, Minovirus, and Salovirus. Lagovirus, Norovirus, Nebovirus, Recovirus, Sapovirus, Valovirus, and Vesivirus infect mammals, while Bavovirus and Nacovirus infect birds, and Minovirus and Salovirus infect fish. Human norovirus is a major cause of acute gastroenteritis in humans. Additionally, unclassified caliciviruses have been found in geese, yellowtail, rainbow trout, lampreys, frogs, and various Australian bird species (14).

In the study, 1,080 viral reads matching caliciviruses were identified. Four nearly complete calicivirus genomes were obtained through software assembly and designated as CALIV1, CALIV4, CALIV5, and CALIV8. The approximate genome lengths of these caliciviruses are 8,856 nt, 6,328 nt, 7,181 nt, and 5,251 nt, respectively, with GC contents of 55.0, 47.9, 55.3, and 47.0%. Partial or complete polyprotein sequences were assembled for these caliciviruses. CALIV1 has a polyprotein length of 2,385 amino acids, CALIV4 has 1,883 amino acids, CALIV5 has 2,127 amino acids, and CALIV8 has 1,524 amino acids (Figure 3A). BLASTp comparisons of the encoded polyproteins on NCBI revealed that CALIV1 and CALIV5 have the highest amino acid sequence identity (89.84 and 90.76%, respectively) with strains isolated from goose fecal samples in China in 2012. While CALIV4 and CALIV8 have the highest amino acid sequence identity (31.9 and 32.68%, respectively) with strains isolated from grey ducks in Australia.

To further determine the genetic evolutionary relationships of the four caliciviruses identified in this study, a phylogenetic tree was constructed based on the polyprotein sequences of CALIV1, CALIV4, CALIV5, CALIV8, and representative strains from different genera within the Caliciviridae family. As shown in Figure 3B, CALIV4 and CALIV8 form a distinct branch, exhibiting a relatively distant genetic relationship from other representative strains. Meanwhile, CALIV1 and CALIV5 cluster together with a virus detected from goose fecal samples in China in 2014 (NC_024078). They do not cluster with representative strains from known genera, suggesting that they may belong to a new genus or species within the Caliciviridae family (Figure 3B).

The Circoviridae family consists of viruses with circular, covalently closed single-stranded DNA (ssDNA) genomes, including the smallest known animal viruses. This family comprises two genera: Circovirus and Cyclovirus. With the development of nested PCR and metagenomic sequencing technologies, the diversity within this family has greatly expanded. Circoviruses have been detected in various human, animal, and environmental samples, indicating their wide distribution in nature. The pathogenicity of many members of this virus family remains poorly understood. Species classification within specific genera is based on genome structure, particularly the position of the replication origin (ori) relative to the coding region. The species demarcation threshold for Circoviridae members is 80% nucleotide sequence identity across the entire genome. Based on this criterion, the Circoviridae family currently includes 101 species (15).

In this study, 15,981 viral sequence reads matching the Circoviridae family were detected from the bar-headed goose sample library. Using the assembly program in Geneious 11.1.2, 10 complete circovirus genomes were obtained and designated as CICV1, CICV3, CICV4, CICV7, CICV10, CICV15, CICV27, CICV36, CICV39, and CICV44. These genomes vary in length, ranging from 3,948 nt to 4,887 nt. Each of these genomes encodes a replication-associated protein. Additionally, CICV15, CICV27, CICV36, and CICV39 encode a capsid protein, while CICV1, CICV3, CICV4, CICV7, and CICV10 encode three putative proteins of unknown function. Notably, CICV44 encodes only a replication-associated protein and lacks a coding frame for the capsid protein, suggesting it is a defective virus (Figure 4A). The GC content of these 10 circoviruses ranges from 29.0 to 51.3%. These newly discovered sequences shared low sequence identity with known sequences in the NCBI database. Specifically, CICV7, CICV36, CICV39, and CICV44 did not find homologous viral sequences. Similar to other circoviruses, the viral replication-associated protein (Rep) contains conserved NTPase/helicase motifs (Walker A, B, and C) (Figure 4B).

Due to the lack of a capsid protein, CICV44 must rely on a helper virus to complete its full lifecycle. When co-infecting host cells, CICV44 competes with the complete virus for the host’s transcription and translation resources. This competition can affect the replication and expression of the complete virus, potentially reducing its infection efficiency. Additionally, the presence of CICV44 can exert evolutionary pressure on the complete virus, selecting for more adaptable viral strains. This selection may alter the antigenic properties and pathogenicity of the virus. By altering the assembly of viral particles or interfering with the host immune response, the defective virus may also help the complete virus to evade the host’s immune surveillance. This evasion can enhance the transmission and persistence of the virus.

To clarify the genetic evolutionary relationships of the 10 circoviruses identified in this study, a phylogenetic tree was constructed using the replication-associated proteins from these strains and representative strains from various genera within the Circoviridae family. As illustrated in Figure 4C, CICV15, CICV27, CICV36, and CICV39 form a distinct branch, clustering with other unclassified circoviruses found in diverse environmental and biological samples. In contrast, CICV1, CICV3, CICV4, CICV7, CICV10, and CICV44 cluster with Po-circo-like viruses. Sequence alignment and phylogenetic analysis confirm that the viruses identified in this study are novel circoviruses.

The Parvoviridae family includes a rich and highly diverse group of viruses. These viruses are characterized by linear ssDNA genomes (4–6 kb) and icosahedral capsids (20–25 nm) (16). The host range includes both invertebrate and vertebrate hosts, and they are primarily divided into three subfamilies: Densovirinae, Hamaparvovirinae, and Parvovirinae. Due to their limited genetic material, most parvoviruses require actively dividing host cells and exhibit host and/or tissue specificity. Some can cause diseases ranging from subclinical to fatal. A few require co-infection with other helper viruses. The Hamaparvovirinae subfamily is a recently established taxon, with members capable of infecting both vertebrates and invertebrates. Chaphamapavirus is a genus within the Hamaparvovirinae subfamily. This virus has been detected in fecal materials from chickens, turkeys, rats, pigs, and bats, as well as in various tissue samples from pigs, including serum, rectal swabs, nasal swabs, and lung lavage fluids. Additionally, a novel avian parvovirus was recently detected in brain, liver, and heart tissues collected from rainbow lorikeets. However, it is currently unclear whether these viruses are associated with any known diseases.

In this study, we identified 18,194 sequence reads matching the Parvoviridae family from the bar-headed goose sample library. Using Geneious 11.1.2 software, we assembled 19 complete genomes of parvoviruses, named PACV2, PACV4, PACV5, PACV6, PACV8, PACV9, PACV10, PACV13, PACV18, PACV22, PACV24, PACV26, PACV28, PACV29, PACV31, PACV34, PACV35, PACV46, and PACV50. These genomes vary in length from 4,348 to 4,705 nucleotides and have GC contents ranging from 37.8 to 40.5%. All identified parvoviruses contain four open reading frames (ORFs) that encode a hypothetical protein, a nucleoprotein, a capsid protein, and a nonstructural protein (Figure 5A). The conserved NTPase/helicase motifs (Walker A, B, B′, and C) were found in the nonstructural proteins of all 19 viruses (Figure 5B). BLASTn analysis in NCBI revealed varying levels of similarity with known chaphamaparvovirus strains. PACV2 and PACV34 share the highest identity with strains isolated from Australian chestnut teal feces in 2018, with 73.84 and 73.77% similarity, respectively. PACV4, PACV6, and PACV18 show similarities of 74.48, 77.54, and 74.60% with strains from mute swan feces in the UK in 2016. PACV5 has 78.08% similarity with a chaphamaparvovirus strain from mute swan feces in the UK in 2016. PACV8 shows 74.87% similarity with a strain from wild bird feces in China in 2016, while PACV9 is 72.87% similar to a strain from black swan feces in China in 2018. PACV10, PACV22, and PACV28 show similarities of 74.06, 73.03, and 75.15% with strains from bar-headed goose feces in China in 2016. PACV13 shares 75.96% similarity with a strain from wild bird feces in China in 2018. PACV24 shows 78.61% similarity with a strain from black swan feces in China in 2018. PACV26 is most similar to a strain from chestnut teal feces in Australia in 2018, with 76.57% similarity. PACV29 shows the highest similarity of 82.32% with a strain from black swan feces in China in 2018. PACV31 and PACV35 are most similar to strains from ruddy shelduck feces in China in 2016, with 79.03 and 73.68% similarity, respectively. PACV46 shows 71.84% similarity with a strain from chicken liver tissue in China in 2022, while PACV50 is 80.59% similar to a strain from mute swan feces in the UK in 2016. These results suggest that the parvoviruses identified in this study are novel, with varying degrees of relatedness to previously known chaphamaparvovirus strains.

To elucidate the genetic evolutionary relationships of the 19 parvovirus strains identified in this study, we constructed a phylogenetic tree based on the NS1 proteins of these parvoviruses and representative strains from the Chaphamaparvovirus genus. The phylogenetic analysis revealed that PACV34, PACV13, and PACV2 clustered together and were closely related to a chaphamaparvovirus strain (NC_077064) isolated from the feces of an Australian chestnut teal in 2018. PACV31 and PACV35 formed an independent branch, while PACV4, PACV5, PACV6, PACV18, and PACV24 clustered with strains isolated from swan feces samples in the UK. PACV9, PACV29, PACV46, and PACV50 also formed an independent branch, showing a distant genetic relationship with other strains, suggesting they might represent new species (Figure 5C). To further analyze the specific species of these newly discovered parvoviruses, we selected multiple reference strains of different species within the Chaphamaparvovirus genus and constructed another phylogenetic tree based on NS1 (Figure 5D). The results indicated that the newly discovered 19 strains do not belong to any existing species, implying they may represent four new species.