Genetic Analysis of Avian Gyrovirus 2 Variant-Related Gyrovirus Detected in Farmed King Ratsnake (Elaphe carinata): The First Report from China

Avian gyrovirus 2 (AGV2), which is similar to chicken infectious anemia virus, is a new member of the genus Gyrovirus. AGV2 has been detected not only in chicken but also in human tissues and feces. This study analyzed 91 samples (8 from liver tissue and 83 from fecal samples) collected from king ratsnakes (Elaphe carinata) from 7 separate farms in Hubei and Henan, China, for AGV2 DNA using PCR. The results demonstrated a low positive rate of AGV2 (6.59%, 6/91) in the snakes, and all the positive samples were collected from the same farm. The AGV2 strain HB2018S1 was sequenced, and its 2376 nt genome comprised three partially overlapping open reading frames: VP1, VP2, and VP3. Phylogenetic analysis revealed that the HB2018S1 and NX1506-1 strains from chickens in China belong to the same clade and that they have a nucleotide identity as high as 99.5%. Additionally, recombination analysis showed that HB2018S1 might originate from the recombination of viruses similar to those detected in chickens and a ferret. A total of 10 amino acid mutation sites (44(R/K), 74(T/A), 256 (C/R), 279(L/Q), and 373(V/A) in AGV2 VP1; 60(I/T), 125(T/I), 213(D/N), and 215(L/S) in AGV2 VP2; and 83(H/Y) in AGV2 VP3) different from those observed in most reference strains were found in the genome of HB2018S1, indicating that the differences may be related to a transboundary movement among hosts, which needs further elucidation. To the best of our knowledge, this study is the first report on an AGV2-infected poikilotherm, suggesting that cross-host transmission of viruses with circular single-stranded DNA genomes would be a public health concern.


Introduction
Avian gyrovirus 2 (AGV2) belongs to the genus Gyrovirus and was originally classified in the family Circoviridae but was reassigned to the family Anelloviridae according to the recent viral taxonomy report of the International Committee on Taxonomy of Virus, which states that AGV2 is not structurally or genetically related to the members of the family Circoviridae [1]. Besides AGV2, the family Anelloviridae includes another recognized species-chicken anemia virus (CAV). CAV was commonly found in chickens worldwide, leading to significant losses in the poultry industry since it was first isolated from the bursa of chicken Fabricius in 1979 in Japan [2]. AGV2 was first detected in diseased chickens in Brazil in 2011, and it harbored approximately 40% nucleotide similarity with that of CAV [3]. The AGV2 genome is approximately 2.38 kb and a circular single-stranded DNA (ssDNA); it mainly comprises three overlapping open reading frames (ORFs) coding for VP1, VP2, and VP3 [4]. AGV2 infections in chickens can result in brain damage, mental retardation, and weight loss [5]. Although other specific symptoms of AGV2 infections have not been confirmed, autopsy-based studies have reported clinical manifestations such as hemorrhage, edema, glandular gastric erosion, and facial and head swelling in infected chickens [5]. Varela et al. (2014) used duplex quantitative real-time PCR assay to assess commercially available poultry vaccines and suggested that the widespread existence of AGV2 is associated with vaccine contamination [6]. AGV2 has also been reported in different regions of Europe, Latin America, Africa, South America, and Asia, indicating its global distribution [7][8][9][10]. In 2011, Sauvage et al. (2011) identified human gyrovirus, which shared 96% nucleotide identity with AGV2, from a skin swab of a healthy individual, indicating that AGV2 may infect humans as well [11]. Furthermore, AGV2 has been also detected in human blood samples [12]. In China, farmed king ratsnake is a highly distributed nontoxic snake species; moreover, snake meat has a high nutritional value and medicinal efficacy. Fortuitously, one farmed king ratsnake that died from bacterial infection was detected as AGV2-positive; therefore, the present study aimed to investigate AGV2 in farmed king ratsnakes in China and amplify its genome. We subsequently performed an in-depth sequence analysis on the basis of genetic evolution and amino acid mutations between the sequenced strain and reference strains.

Samples and Virus Detection
In 2018, 91 samples (8 from liver tissues of snakes that died from bacterial infection and 83 from feces collected with individual sterile swabs and placed in sterile collection tubes) collected from 7 king ratsnakes from different farms (farms A-G) in Hubei and Henan, China, were tested for AGV2 by PCR using AGV2-specific primers (AGV2-F (5 -CGTGTCCGCCAGCAGAAAC-3 ) and AGV2-R (5 -GGTAGAAGCCAAAGCGTCCAC-3 ); from nucleotides 656-1001) on the basis of the highly conserved region of the AGV2 genome, as described previously [12]. The liver tissues were obtained via dissection to avoid contamination, and the liver tissue and fecal samples were washed with 1 mL phosphate-buffered saline, frozen and thawed thrice, immersed in liquid nitrogen, and ground into a homogenate. The homogenate was centrifuged at 5000 rpm for 10 min, and 0.2 mL of the supernatant was used for nucleic acid extraction. DNA and RNA were extracted using a DNA/RNA extraction kit (TransGen Biotech, Beijing, China) according to the manufacturer's instructions. The extracted DNA was quantified using the NanoDrop ND-1000 (Thermo Scientific, Waltham, MA, USA) according to the manufacturer's instructions using 1 µL of the DNA sample, and the quantified extracted DNA (concentration: >100 ng; ratios of 260:280: >1.8, and ratios of 260:230: >2.0) was stored at −20 • C until use.

Recombination and Phylogenetic Analysis
SeqMan (DNASTAR, Lasergene, Madison, Wisconsin) was used for contig-assembly of the partial sequences. The whole-genome sequence of the snake-originated strain was submitted to GenBank under the accession number MK840982. After alignment of the HB2018S1 and 25 AGV2 reference sequences downloaded from GenBank (detailed information of each reference strain is presented in Supplementary Materials Table S1) using Clustalx v1.83, phylogenetic trees were constructed on the basis of the nucleic acids of whole genomes using MEGA v.6.0 with the neighbor-joining method using the Kimura 2-parameter model, gap pairwise deletion, and bootstrap analysis of 1000 replicates [13]. Furthermore, the RDP 4.36 and SimPlot 3.5.1 software were used for the genetic recombination analysis of AGV2, and the assumed parental strain was used as the reference strain. RDP was executed with default parameter settings, which include the RDP, GENECONV, and Maxchi. The results obtained by RDP4 were confirmed using a boot-scanning method in the SimPlot program v.3.5.1 [14]. Meanwhile, phylogenetic trees were also constructed as described above on the basis of the genome region relative to the predicted recombination events. In addition, mutation sites were compared according to the amino acid sequences of the three ORFs.

Positive Rates of AGV2 in the Snakes in China
To increase our understanding of AGV2 infection in farmed king ratsnakes in China, 91 samples were collected and analyzed. The results indicated that five fecal samples (6.02%, 5/83) and only one liver tissue sample (12.5%, 1/8) were positive for AGV2. All six of these positive samples were collected from the same snake farm in Hubei, China (the detection results of AGV2 in this farm are presented in Supplementary Materials Figure S1). Detailed information on the snake samples positive for AGV2 is listed in Table 1.

Whole-Genome Sequencing of AGV2 in Snakes
On the basis of the sequencing and contig-assembly results, the AGV2 strain identified in the liver tissues of king ratsnakes was named HB2018S1, and its genome comprised a 2376 nt long sequence. The sequencing results showed that the genomes of the six AGV2-positive snakes were 100% identical, indicating that they were the same strain. The nucleotide identity of HB2018S1 with the same clusters of the NX1506-1, NX1506-2, NX1510, and JL1508 reference strains was up to 99.5%. The lowest nucleotide identity of 95.3% was detected between HB2018S1 and S53/It from chickens. In addition, the nucleotide similarities between HB2018S1 and the three human strains 915F06007FD, CL33, and JQ690763 were 95.6%, 97.1%, and 96.8%, respectively.

Phylogenetic Analysis
Phylogenetic analysis revealed that the AGV2 and 25 full-length genome sequences downloaded from GenBank could generally be sub-grouped as two branches ( Figure 1). HB2018S1 and most of the downloaded sequences belonged to branch I, which was subsequently divided into three clusters and marked.

Phylogenetic Analysis
Phylogenetic analysis revealed that the AGV2 and 25 full-length genome sequences downloaded from GenBank could generally be sub-grouped as two branches ( Figure 1). HB2018S1 and most of the downloaded sequences belonged to branch I, which was subsequently divided into three clusters and marked. The SimPlot 3.5.1 and RDP 4.36 software were used to predict the possibility of recombination events in the AGV2 genome. As shown in Table 2, four recombination events were identified. In recombination event 1, the recombination sites of seven strains were identical, and all the strains covered belonged to branch I; moreover, in this study, HB2018S1 was also included in recombination event 1. The recombination occurred between the relatives of viruses within the clade of G13 (from CL33 to 915F06007FD) and the relatives of viruses within the clade of HLJ1603-2 (from Ave3 to HLJ1603-1). The recombinant virus gave origin to the whole clade that includes HB2018S1 ( Figure 2). The recombination event was mapped at 100-1542 nts, and the region related to recombinant event is shown in Figure 3, which was created using the SnapGene software. This region (100-1542 nts) was similar to G13 (99.5%), and the remaining AGV2 genes were similar to HLJ1603-2 (99%). To further The SimPlot 3.5.1 and RDP 4.36 software were used to predict the possibility of recombination events in the AGV2 genome. As shown in Table 2, four recombination events were identified. In recombination event 1, the recombination sites of seven strains were identical, and all the strains covered belonged to branch I; moreover, in this study, HB2018S1 was also included in recombination event 1. The recombination occurred between the relatives of viruses within the clade of G13 (from CL33 to 915F06007FD) and the relatives of viruses within the clade of HLJ1603-2 (from Ave3 to HLJ1603-1). The recombinant virus gave origin to the whole clade that includes HB2018S1 (Figure 2). The recombination event was mapped at 100-1542 nts, and the region related to recombinant event is shown in Figure 3, which was created using the SnapGene software. This region (100-1542 nts) was similar to G13 (99.5%), and the remaining AGV2 genes were similar to HLJ1603-2 (99%). To further prove the occurrence of recombination events, we split the genome into two parts at recombinant sites 100 and 1542 nts and constructed phylogenetic trees. As shown in Figure 4A, according to the 100-1542 nts of the whole genome, HB2018S1 and G13 belong to the same branch and are far from HLJ1603-2. Conversely, phylogenetic trees constructed on the basis of the remaining partial whole genome show that HB2018S1 and HLJ1603-2 belong to the same branch and are far from G13 ( Figure 4B). These results were in accordance with the recombination prediction.

Analysis of the Mutated Amino Acids
The AGV2 genome comprises three overlapping ORFs, which encode VP1, VP2, and VP3 proteins. ORF3, which is located at sites 953-2335 and encodes VP1, comprises 460 amino acid residues. Three replication motifs were observed in the VP1 of HB2018S1: FAALS (325-329 residues), RRWLTLV (363-369 residues), and KAMA (412-415 residues). VP1 of HB2018S1 was 97.4-98.9% identical with the reference strains in terms of deduced amino acids. On the basis of the sequence alignment between HB2018S1 and the reference strains, 20 amino acid mutations were observed for VP1 (displayed in gray in Table 3 The ORF1 encoding AGV2 VP2 was located at sites 450-1145 and encoded a 231-amino acid long protein sequence, and the VP2 sequences of HB2018S1 were 94.0%-98.3% identical to those of the reference strains in terms of deduced amino acids. Sequence alignment indicated 15 hyper-mutated amino acid sites located in the AGV2 VP2 of HB2018S1 and reference sequences, whereas G17 from ferrets had an S insertion at site 162. In addition, the substitutions 60I/T, 125T/I, 213D/N, and 215L/S occurred only in HB2018S1 (displayed in gray in Table 4).
a Mutations were displayed in gray.

Substitution of Amino Acid
The untranslated region (UTR) of HB2018S1 was located between the canonical polyadenylation site (AATAAA) and the start site of transcription, which included six closely matching direct repetition (DR) regions with a length of 22 nts. Among these, five DRs were identical (5 -GTACAGGGGGGTACGTCACCAT-3 ). The different DR was the last DR. The three bases at its 3 end were "A, G, and C" of HB2018S1, which were similar to the reference strains.

Discussion
To the best of our knowledge, AGV2 has been detected in three different hosts-chickens, humans, and ferrets [10][11][12]. Most AGV2 strains are mainly detected in chickens, which is the original host. In the present study, AGV2 was detected, for the first time, in cold-blooded farmed king ratsnakes; however, how the snakes were infected with the virus remains unknown. Farmed snakes are usually given a specific feed and occasionally raw meat. In this study, snakes from only one farm were positive for AGV2. According to a survey regarding the feed, the infected snakes (all from the same farm) had been fed chicken meat 1 month before sampling, whereas those in the other farms had been fed chicken meat supplied from different chicken farms within or after 1 week. On the basis of a random investigation of all the poultry farms that supplied chicken meat, the positive ratio for AGV2 of each farm reached 40%. It is supposed that the AGV2 detected in snakes originated from infected chickens, and that HB2018S1 might harbor the capability of cross-host transmission for only one farm whose snakes were tested to be AGV2-positive. Although it is challenging to prove that the origin of the AGV2 detected in the snakes was the chicken meat used as feed because the virus was not successfully isolated from chicken meat, AGV2 detected in snake liver tissues may actually indicate that the virus replicates in the host. Furthermore, the positive liver tissue and fecal samples collected from the same farm suggested that the fecal samples positive for AGV2 did not correspond to residual viruses from the gut because the snakes from other farms negative for AGV2 may also have been fed infected chickens, given that the AGV2 infection rate in chickens is very high in China [5,12]. Whether snakes are infected with AGV2 because of alimentary transmission following chicken consumption requires further investigation. In 2012, AGV2 was detected in the feather sacs of two adjacent chicken flocks in southwestern Brazil, suggesting a horizontal transmission mode of AGV2; however, more studies are required to demonstrate this phenomenon. From 2011 to 2013, some researchers conducted an epidemiological analysis of AGV2 in the serum and feces of some individuals to determine whether AGV2 could be replicated in humans. Whether the human AGV2 is attributable to metabolic residues in poultry food remains to be clarified. Phylogenetic analysis showed that the genome sequence of HB2018S1 was highly similar to those of NX1506-1, NX1506-2, and NX1510 from China. However, the local regions, wherein these strains were identified, are far from the collection site of HB2018S1, and the genome sequences of some strains from the chicken farms that supplied the chicken meat were distant from HB2018S1 (data not shown). Therefore, the AGV2 transmission mode requires further investigation.
The AGV2 genome is 2376 nt long, similar to the AGV2 genome from chickens in China, having three overlapping ORFs, similar to those in chicken infectious anemia virus (CIAV) encoding VP1, VP2, and VP3. At the protein/amino acid level, the similarity of the VP1, VP2, and VP3 proteins in AGV2 to the VP1, VP2, and VP3 proteins in CIAV was 38.8%, 40.3%, and 32.2%, respectively [3]. The UTR, which is located between the polyadenylation site and the site of transcription initiation, is similar to CIAV in structure [15,16].
Previous reports have indicated that the AGV2 infection rate in chickens increases with time [5,12]. If AGV2 was transmitted from chickens to snakes in this study, the most probable explanation for only one virus strain or virus originating from one strain being detected and sequenced in snakes was the mutation of functional genes, resulting in trans-host infection. The VP1 protein of CIAV is the only capsid protein contained in the virus, and constitutes the antigen that elicits the appearance of neutralizing antibodies. There are three replication motifs on the VP1 sequence in CIAV, including FATLT (at 313-317 residues), QRWHTLV (at 351-357 residues), and YALK (at 402-405 residues) located on the VP1 sequence conserved in CIAV [17]. Similar sequences were observed in the VP1 of HB2018S1, and the sites are completely conserved in the currently reported AGV2 sequences. On the basis of part of the AGV2 sequences, Yao et al. (2016) hypothesized that a high-variant region exists in the middle and back of the VP1 protein, ranging from 288 to 314 residues. In the present study, only eight strains of viruses, including GS1512, HLJ1508, G13, JL1511, HM590588, JQ690763, RS/BR/15/25, and G17, had mutations, and the mutation sites were not identical. Because of the limited number of sequences, it was not possible to predict the high-variation regions accurately. However, further studies are required to verify whether the five mutations in the VP1 protein in HB2018S1, which are different from those in the other sequences, may affect viral host selectivity.
The VP2 protein in CIAV is a dual-specific protein phosphatase [18]. The WX7HX3CXCX5H sequence of the phosphatase motif is also highly conserved in CAV, torque teno virus (TTV), and TTV-like mini virus (TLMV) [17,19]. A similar motif in the VP2 protein of AGV2 are WLRQCARSHDEICTCGRWRSH (95-115) residues. Site-directed mutagenesis of VP2 revealed that VP2 affects the replication of viral particles, resulting in cytopathological changes and downregulating the levels of major histocompatibility complex I (MHC I) in infected cells [18]. VP1 and VP2 are necessary for viral particle replication. Noteborn [20] also observed interactions between VP1 and VP2 proteins on the basis of an immunoprecipitation test, which further confirmed that the nonstructural protein VP2 plays a scaffold role in virus particle assembly. In the VP2 protein, HB2018S1 has four different amino acid mutations and an S insertion site of G17 from ferrets. Further studies are warranted to determine whether such mutations alter VP2 function and, in turn, affect virus replication, and whether they affect the scaffold function of the VP2 protein.
The VP3 protein in CIAV is also known as apoptin. In normal cells, the VP3 protein exists in the cytoplasm. However, in cancerous cells, it exists in the nucleus and induces apoptosis [21,22]. It has also been demonstrated that the VP3 protein in AGV2 induces specific apoptotic functions in tumor cells, whereas the change in molecular function caused by a mutation at site 83H/Y of the VP3 protein in HB2018S1 still requires further experimental evidence.
In the present study, in addition to the AGV2 identified from chickens, humans, and ferrets, an AGV2 strain from a cold-blooded farmed snake was identified, which could provide insights into the modes of transmission of AGV2 and its host diversity.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-0817/8/4/185/s1, Table S1: The sequence information of reference strains used in this study, Figure S1: Detection results of AGV2 in Farm B samples by electrophoresis.
Author Contributions: J.J., X.X. and Q.X. designed the study and wrote the manuscript; Q.W. and Q.C. performed the sampling and data collection; Y.K. performed the clinical investigations, necropsies, and molecular tests; L.Y. performed the molecular genetic studies and helped to draft the manuscript. Q.W. and X.X. contributed equally to this work All the authors have read and approved the final manuscript.