What is for dinner? Viral metagenomics of US store bought beef, pork, and chicken

We describe here the metagenomics-derived viral sequences detected in beef, pork, and chicken purchased from stores in San Francisco. In beef we detected four previously reported viruses (two parvoviruses belonging to different genera, an anellovirus, and one circovirus-like virus) and one novel bovine polyomavirus species (BPyV2-SF) whose closest relatives infect primates. Detection of porcine hokovirus in beef indicated that this parvovirus can infect both ungulate species. In pork we detected four known parvoviruses from three genera, an anellovirus, and pig circovirus 2. Chicken meat contained numerous gyrovirus sequences including those of chicken anemia virus and of a novel gyrovirus species (GyV7-SF). Our results provide an initial characterization of some of the viruses commonly found in US store-bought meats which included a diverse group of parvoviruses and viral families with small circular DNA genomes. Whether any of these viruses can infect humans will require testing human sera for specific antibodies.


Introduction
Viral metagenomics characterizes known and divergent viruses in sample preparations enriched for viral nucleic acids. Host DNA or RNA concentration is reduced by nuclease treatment while viral nucleic acids are protected within filterable viral particles (Allander et al., 2001). Following non-specific amplification and deep sequencing the viral sequences are recognized through their translated protein sequence similarity to known viral proteins. Here we characterize viral sequences in beef, pork, and chicken bought from stores in San Francisco. Such viruses are of interest given the enormous scale of their exposure to humans consuming meat products.
Meat and meat products are the source of numerous enteric infections in humans. Enteric viruses harbored by farmed animals may contaminate meat samples during slaughter and meat preparation (Machnowska et al., 2014). Hepatitis E virus (HEV) replicating in domesticated or wild-pigs can also be transmitted to human through undercooked pork consumption (Christou and Kosmidou, 2013;Meng, 2013). At slaughter time (22-29 weeks), HEV was present in 12% and 41% of the plasmas and stools from Canadian pigs, respectively (Leblanc et al., 2007). Viral transmission between farm animal species or between domesticated and wild animals may also occur increasing risk of adaptation to new host and changes in viral tropism (Ma et al., 2009;Sonnberg et al., 2013). Transmission of viruses in meat or blood may also occur through non-enteric routes such as skin cuts occurring during hunting and butchering as seen for simian T-lymphotropic and simian foamy virus (Betsem et al., 2011;Gessain et al., 2013;Switzer et al., 2004). Cross-species transmissions, likely due to hunting and/or butchering non-human primates, initiated various HIV epidemics (HIV-1 groups M, N, O, P and HIV-2 groups A-H) (Sharp and Hahn, 2011) and Ebola virus outbreaks (Leroy et al., 2004). The possibility of bovine viruses playing a role in carcinogenesis has been raised based on the association of colorectal cancer with consumption of under-cooked beef (Zur Hausen, 2012).

Overview of virome
The Illumina MiSeq 2 Â 250 bases run generated a total of 4,944,932 unique sequence reads. Sequence reads were assembled de novo and compared to the GenBank non-redundant protein database using BLASTx. 475 sequence reads showed similarity to known eukaryotic viral sequences with an E-value cutoff 10 À 5 .  Table 1 presents the characterization of the eukaryotic viral sequences detected in the meat samples.
In the beef samples we detected eukaryotic viral sequences related to parvoviruses, polyomaviruses, anelloviruses and circoviruses (Table 1). In the pork samples we identified viral sequences related to parvoviruses, anelloviruses and circoviruses (Table 1). In chicken diverse gyroviruses were detected (Table 1).

Parvoviruses in beef samples
Sequences identical to the previously described bovine parvovirus 3 (BPV3: prototype of the ungulate erythroparvovirus 1 species) were identified. Here we used the recently described taxonomy for the family Parvoviridae (Cotmore et al., 2014). BPV3 was one of two parvoviruses previously identified from bovine serum (Allander et al., 2001). Sequences identical to porcine hokovirus (P-PARV4 in the ungulate tetraparvovirus 2 species), first sequenced from porcine tissues (Lau et al., 2008), were also identified in beef. Porcine hokovirus is distinct from bovine hokovirus (B-PARV4-1 in the distinct ungulate tetraparvovirus 1 species) with only 63-65% nucleotide identity between them in ORF1 and ORF2 (Lau et al., 2008).
In order to confirm the existence of porcine hokovirus in beef samples, PCR screening was performed using two sets of nested primers. Twenty-two beef samples, including the 10 samples analyzed here and 12 beef samples collected from stores in San Francisco in 2010 , were tested. PCR results indicated two beef samples were positive, one from a pool of five samples and the other from a beef sample collected in 2010. The PCR products were sequenced confirmed, sharing 99% sequence identity with a porcine hokovirus strain in GenBank (JQ700068).

A novel bovine polyomavirus
Polyomaviruses are small, circular DNA viruses that can cause persistent infections in numerous mammals and can be oncogenic in some hosts. The one previously reported bovine polyomavirus (BPyV) was originally isolated from bovine serum (Schuurman et al., 1990). Due to its high prevalence in cows BPyV is also a frequent environmental contaminant (Hundesa et al., 2006). A high seroprevalence to BPyV has been reported in people with high exposure to cattle (Parry and Gardner, 1986) but no specific illness in either cows or humans has been attributed to BPyV.
In one beef samples pool 31 sequence reads generating eight different contigs were found to have strong amino acid similarity to polyomaviruses. PCR screening of the individual samples with a set of nested primers designed on the LT-Ag sequence indicated that a single beef sample was positive. Using inverse PCR the complete genome was acquired and Sanger sequenced (GenBank KM111535). The sequences of the primers used for PCR screening and inverse PCR are shown in "Materials and methods".
The circular genome of the bovine polyomavirus (named BPyV2-SF) was 5085 bases encoding five major proteins including an early region with regulatory ST-Ag and LT-Ag and a late region coding for the VP1, VP2, and VP3 structural proteins (Fig. 1A), and a 214 aa alternative Large T ORF (ALTO) protein lying in the second exon of LT-Ag gene. The early region and late region were separated by a regulatory region that contained typical polyomavirus features (Pipas, 1992) including an AT-rich region on the left side of the putative replication origin. Four repeats of the consensus pentanucleotide LT-Ag binding site GAGGC were present. In polyomaviruses VP3 is translated by a leaky ribosomal scanning starting at an internal methionine codon of the larger VP2 ORF. The known polyomaviruses can be divided into "VP3 þ " species with a clear homolog of the VP3 N-terminal MALXXΦ motif (Φ is an aromatic residue) and "VP3-less" species lacking the motif (Schowalter and Buck, 2013). Phylogenetic analysis based on complete genome of known polyomaviruses revealed that VP3less species cluster together into a monophyletic clade to which BpyV2-SF also belongs (Fig. 1B). Sequence analysis reveals that BpyV2-SF fits the characteristics of VP3-less polyomavirus species, such as having an ALTO protein within the LT-Ag frame, lacking the VP3 N-terminal MALXXΦ motif, and missing $80 aa near the carboxyl terminus of the VP2 compared with the longer VP3 þ species (Schowalter and Buck, 2013).
To determine the divergence in sequence between BPyV2-SF and other VP3-less polyomavirus species, amino acid sequence alignments of LT-Ag, VP1, and VP2 were performed, and neighborjoining trees generated. The sequence alignments include BPyV2-SF, 10 representative VP3-less polyomavirus strains, and seven representative VP3 þ strains. Phylogenetic analysis indicated that BpyV2-SF is highly divergent from these other polyomaviruses and that it fell within the VP3-less polyomavirus monophyletic clade together with polyomaviruses from gorilla, chimpanzee, humans, eidolon (megabat), and raccoon (Fig. 1B). Table 2 presents the amino acid distance between BPyV2-SF and the closest relative in LT-Ag, ST-Ag, VP1 and VP2. The high divergence of protein sequences confirms that BPyV2-SF is a novel polyomavirus species.

Circovirus-like sequence in the beef
In a beef pools two sequence reads showed high sequence identity with a previously reported circovirus-like virus PorkNW2 sequence (GenBank accession HQ738638), previously found in both pork and beef samples from San Francisco stores using Circoviridae consensus primers . PCR screening with a set of  Another ORF potentially encoding 243 aa protein arranged in the same coding direction as Rep protein was also found yielding an ORF map similar to that of NW2 ( Fig. 2A and B). BLASTp search indicated part of the putative protein (60 aa) shared 47% amino acid sequence identity with the M169 protein of unknown function of Murid herpesvirus (GenBank no. HE610454).

Viruses in pork
We identified sequences from pork-derived viral particles that were identical to porcine parvovirus 4 (PPV4 in the ungulate copiparvovirus 2 species). PPV4 was initially identified in a lung lavage from a pig co-infected with pig circovirus 2 (Cheung et al., 2010). We also detected ungulate tetraparvovirus 2 (aka porcine hokovirus, the same virus detected in beef samples) initially reported to infect pigs in Honk Kong at high frequency (Lau et al., 2008). Sequences from a third parvovirus species in a third genus (ungulate protoparvovirus 1) were also detected. Last, a fourth set of parvovirus-like sequences closest to a still unclassified copiparvoviruses labeled porcine parvovirus 6 (PPV6 deposited March 31, 2014 GenBank NC_023860) from pigs in China were also detected. Because PPV6 NS1 shows closest identity to PPV4 but with only 58% identity it likely represents a new copriparvovirus species (proposed name ungulate copiparvovirus 3) (Cotmore et al., 2014). Four parvovirus species in three distinct Parvoviridae genera were therefore detected in these few pork samples.
Also detected in both pork and beef samples were anelloviruses which have been reported in the blood of numerous mammals including cows and pigs (Brassard et al., 2008). Pig circovirus 2 (PCV2) a frequent infection of pigs occurring world-wide was present.

Gyroviruses in chicken
Gyroviruses are small nonenveloped DNA viruses with a negative sense single-stranded circular DNA of about 2.3 kb. The genus Gyrovirus is currently in the Circoviridae family but closer (although still distinct) in genetic organization to members of the Anelloviridae (Hino and Prasetyo, 2009). Here sequences from six different gyroviruses were detected in the chicken meat samples: chicken anemia virus (CAV), human gyrovirus 1 (HGyV1) and the very closely related gyrovirus 2 (GyV2), GyV3, GyV4 and a novel gyrovirus here named GyV7. BLASTx search indicated that the known gyroviruses found here shared high protein sequence identities (492%) with their best matches in the GenBank, while the novel gyrovirus (GyV7-SF) showed only 57% protein sequence similarity to its closest matches (Table 1). PCR screening with a set of nested primers designed on the novel sequences indicated that only one chicken sample was positive for GyV7-SF. The complete circular genome of GyV7-SF1 was acquired using inverse PCR and Sanger sequencing (GenBank KM111536). The GyV7-SF genome is 2439 bases, with a 514 base non-translated region (NTR) and a polyadenylation signal [AATAAA]. ORF prediction revealed three overlapping ORFs for VP1 (465-aa), VP2 (238-aa) and VP3/Apoptin (130-aa). The genome organization of GyV7-SF was similar to that of the other gyroviruses (Fig. 3A). BLASTp search based on VP1, VP2, and VP3 all revealed the best match with GyV6 (NC_022789), with protein sequence identities of 49%, 53%, and 63%, respectively. Phylogenetic based the three major protein sequence all indicated that GyV7-SF clustered loosely with GyV6 ( Fig. 3B-D), which was detected in diarrhea samples from Tunisian children (Gia Phan et al., 2013).

Discussion
Our results indicate that diverse parvoviruses are readily found in meat samples from hoofed farm animals. Parvoviruses are frequently detected in domesticated and wild animal as well as human feces (Bodewes et al., 2013;Canuti et al., 2011;Handley et al., 2012;Sharp et al., 2010;Väisänen et al., 2014). Some parvovirus can be highly pathogenic, notoriously canine parvovirus which recently adapted from cats to dogs possibly through intermediate species (Stucker et al., 2012) and is thought to be readily transmitted across carnivore species (Allison et al., 2013). Parvoviruses can also be detected in blood as exemplified by PARV4 (primate tetraparvovirus 1) and parvovirus B19 (primate erythroparvovirus 1) which is typically a benign infection but which can cause serious anemia or fetal hydrops in susceptible individuals or when infecting fetus (Feldman et al., 2010;Morinet et al., 2011). Parvovirus B19 can also be deposited in muscles and other tissues where it can be PCR amplified even decades after primary infection (Norja et al., 2006). The detection of ungulate erythroparvovirus 1 virus in beef may therefore reflect either active replication or particle deposition into tissues. Porcine hokovirus (ungulate tetraparvovirus 2) was also found in beef indicating that the tropism of this viral species may include both pigs and cattle.
As expected from its previously reported high frequency in pigs (Lau et al., 2008) porcine hokovirus was also found in the pork samples analyzed here. The lower number of sequence reads in pork than in beef (Table 2) indicates that laboratory contamination from the pork to the beef samples was unlikely. The possibility that the deep sequencing and PCR detection of parvoviruses in beef or pork were due to cross-contamination during slaughter, shipping, or on store shelves cannot be excluded and will require confirmation in independent studies.
Pork samples also contained parvovirus sequences from three genera which had been previously described in other pig samples. Surprisingly we did not detect any sequences belonging to the Bocaparvovirus genus frequently detected in pig fecal samples (Jiang et al., 2014;Shan, Lan, et al., 2011) possibly due to a tropism generally limited to respiratory or gut tissues.
The presence of porcine circovirus 2 in pork may be due to either inoculation with a live attenuated vaccine frequently administered to piglets or from a natural infection (Zhai et al., 2014). A recent study did not detect human antibodies to PCV2 despite frequent consumption of infected pork therefore reducing concern about zoonosis to humans (Burbelo et al., 2013;Victoria et al., 2010). The 0.8 kb circular genome found in beef reassembled a 1.2 kb genome previously sequenced in both pork and beef from the US . Short, PCV2-like, rearranged genomes (483-896 bases) have also been reported in Chinese pigs Wen, He, Ni, et al., 2012;Zhai et al., 2014). The generation of such truncated and rearranged genomes may be related to that seen when viruses are passaged at high multiplicity of infection leading to defective interfering viral particles retaining only essential Cis elements required for genome replication and requiring transcomplementation from a helper virus (Mankertz et al., 1997). A related phenomenon may be the generation of truncated ssDNA circular geminivirus and anellovirus genome when passaged in vitro (De Villiers et al., 2011).
Chicken meat showed the presence of numerous gyroviruses including the highly prevalent CAV reported world-wide (Schat, 2009;Toro et al., 2006). The gyrovirus genus includes several species, namely avian gyrovirus 2 (AGV2) from chicken (Rijsewijk et al., 2011), the closely related HGyV1 from human skin (Sauvage et al., 2011), GyV3 and GyV4 from children feces and chicken meat (Chu et al., 2012;Phan, Li, et al., 2012), and GyV5 and GyV6 from children feces (Gia Phan et al., 2013), many of which were also detected in chicken and children feces from South Africa (Smuts, 2014). The characterization of a seventh species of gyrovirus (GyV7) in this limited sampling of chicken indicates that still more gyroviruses from chicken and other birds are likely to be identified. The current lack of report describing gyrovirus seroconversion in human or replication in mammalian cells indicates that a possible explanation for the frequent detection of gyrovirus DNA in human feces (gyrovirus species GyV3-GyV6) (Chu et al., 2012;Gia Phan et al., 2013;Phan, Li, et al., 2012) and recently cat feces (CAT-CAV) (X.  may reflect their dietary origin from eating chicken or other birds. Only DNA viral sequences were identified in this limited sampling of muscle tissues from presumably healthy animals while RNA viruses such as picornaviruses, caliciviruses, coronaviruses, and astroviruses tend to predominate in the feces of farm animals (Sachsenröder et al., 2014;. The absence of detected RNA viruses in the muscle tissues tested here may reflect a low level of contamination with fecal material, a general preference of RNA viruses for enteric rather than systemic infections, and/or as recently reported (Sachsenröder et al., 2014) the older age of the animals used for meat production.

Sample preparation and metagenomic analysis
Raw pork (7 pieces and 3 ground), beef (7 pieces and 3 ground), chicken (6 pieces and 4 ground) were purchased from six randomly selected meat markets in the Sunset district of San Francisco, USA from Jan 12-28, 2014. Two samples of beef, pork, and chicken were purchased from each of five stores except for two stores in which only a single sample of chicken was purchased. The remaining two chicken samples were purchased from a sixth store. All ground meats were purchased already wrapped in plastic labeled with name of national conglomerates. Neat piece of meats was collected and immediately placed in zippered plastic bags. All samples were frozen the same day at À 80 1C until further processing. None of the meat was butchered or ground in the shops. In order to avoid potential contamination, small inner section of the meat pieces was cut out and used. Meat pieces ($25 mg) were frozen and thawed on dry ice three times, homogenized, filtered, and nuclease treated as previously described to enrich for nucleic acids within viral particles . Remaining total nucleic acid was then isolated using MagMAX™ Viral RNA Isolation Kits (Qiagen) according to manufacturer's protocol.
A library of randomly amplified PCR products from viral RNA and DNA was prepared by using a Nextera DNA Sample Prep Kit and sequenced on the MiSeq Illumina platform (Illumina, San Diego, CA, USA). Six separate pools of nucleic acids from the 30 meat specimens, each pool contained five specimens from the same farm animal species, were generated. The viral nucleic acid pools, containing both DNA and RNA viral sequences, were subjected to reverse transcription reactions with SuperScript III reverse transcriptase (Invitrogen) and 100 picomoles of random hexamer primer, followed by a single round of DNA synthesis using Klenow fragment polymerase (New England BioLabs). A library was then constructed using Nextera XT DNA Sample Preparation Kit (Illumina) and then sequenced using the Miseq Illumina platform with 250 bases paired ends with a distinct molecular tag for each pool. Resulting raw reads were trimmed for quality and primer, and de novo assembled into contigs. Sequences and contigs were compared to the GenBank non-redundant protein database using BLASTx with an E-value cutoff of o10 À 5 .
PCR screening and genome sequencing PCR screening was performed for the bovine polyomavirus, circovirus-like genome, the porcine hokovirus in beef, and the novel gyrovirus in chicken. Inverse PCR was used to generate the complete genome of the novel polyomavirus and the circoviruslike sequence in beef. PCR to bridge sequence gaps was used to acquire the whole genome of the novel gyrovirus in chicken. Sequences and characteristics of the primers used in PCR screening and inverse PCR are shown in Supplemental Table 1. Sanger method was used for sequencing complete viral genomes and to confirm the sequence of PCR products.

Phylogenetic analysis
Phylogenetic analyses were performed based on the predicted amino acid or nucleotide sequences in the present study, their closest viral relatives based on best BLASTx hits, and representative members of related viral species or genera. Sequence alignment was performed using CLUSTAL W with the default settings. Aligned sequences were trimmed to match the genomic regions of the viral sequences obtained in the study. A phylogenetic tree with 100 bootstrap resamples of the alignment data sets was generated using the neighbor-joining method based on the Jones-Taylor-Thornton matrix-based model in MEGA5.0 (Tamura et al., 2011). Bootstrap values (based on 100 replicates) for each node are given. Putative ORFs in the genome were predicted by NCBI ORF finder. Putative exon and intron were predicted by Netgenes2 at http:// www.cbs.dtu.dk/services/NetGene2/.

Nucleotide sequence accession numbers
The genomes of known and new viruses described in detail here were deposited in GenBank under the following accession numbers: KM111535-KM111537. The raw sequence reads from the metagenomic library were deposited in the Short Read Archive of GenBank database under accession number SRX652341.