PLoS OnePLoS ONEplosplosonePLoS ONE1932-6203Public Library of ScienceSan Francisco, USA229371423427170PONE-D-12-0803510.1371/journal.pone.0043986Research ArticleBiologyComputational BiologyGenomicsMetagenomicsGenomicsGenome SequencingMetagenomicsMicrobiologyVirologyZoologyMammalogyVeterinary ScienceVeterinary DiseasesVeterinary VirologyIdentification of a Novel Bat Papillomavirus by MetagenomicsMetagenomics for a Novel Papillomavirus DiscoveryTseHermanTsangAlan K. L.TsoiHoi-WahLeungAndy S. P.HoChi-ChunLauSusanna K. P.WooPatrick C. Y.YuenKwok-Yung*State Key Laboratory of Emerging Infectious Diseases, Department of Microbiology, The University of Hong Kong, Queen Mary Hospital, Pok Fu Lam, Hong Kong Island, Hong KongDavisToddEditorCenters for Disease Control and Prevention, United States of America* E-mail: kyyuen@hkucc.hku.hk

Competing Interests: The authors have declared that no competing interests exist.

Conceived and designed the experiments: HT KYY. Performed the experiments: AKLT HWT CCH ASPL. Analyzed the data: SKL PCYW. Wrote the paper: HT KYY.

2012248201278e4398683201227720122012Tse et alThis is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

The discovery of novel viruses in animals expands our knowledge of viral diversity and potentially emerging zoonoses. High-throughput sequencing (HTS) technology gives millions or even billions of sequence reads per run, allowing a comprehensive survey of the genetic content within a sample without prior nucleic acid amplification. In this study, we screened 156 rectal swab samples from apparently healthy bats (n = 96), pigs (n = 9), cattles (n = 9), stray dogs (n = 11), stray cats (n = 11) and monkeys (n = 20) using a HTS metagenomics approach. The complete genome of a novel papillomavirus (PV), Miniopterus schreibersii papillomavirus type 1 (MscPV1), with L1 of 60% nucleotide identity to Canine papillomavirus (CPV6), was identified in a specimen from a Common Bent-wing Bat (M. schreibersii). It is about 7.5kb in length, with a G+C content of 45.8% and a genomic organization similar to that of other PVs. Despite the higher nucleotide identity between the genomes of MscPV1 and CPV6, maximum-likelihood phylogenetic analysis of the L1 gene sequence showed that MscPV1 and Erethizon dorsatum papillomavirus (EdPV1) are most closely related. Estimated divergence time of MscPV1 from the EdPV1/MscPV1 common ancestor was approximately 60.2–91.9 millions of years ago, inferred under strict clocks using the L1 and E1 genes. The estimates were limited by the lack of reliable calibration points from co-divergence because of possible host shifts. As the nucleotide sequence of this virus only showed limited similarity with that of related animal PVs, the conventional approach of PCR using consensus primers would be unlikely to have detected the novel virus in the sample. Unlike the first bat papillomavirus RaPV1, MscPV1 was found in an asymptomatic bat with no apparent mucosal or skin lesions whereas RaPV1 was detected in the basosquamous carcinoma of a fruit bat Rousettus aegyptiacus. We propose MscPV1 as the first member of the novel Dyolambda-papillomavirus genus.

This work is partly supported by the Research Grant Council, University Development Fund and Strategic Research Theme Fund, The University of Hong Kong; the Tung Wah Group of Hospitals Fund for Research in Infectious Diseases; the HKSAR Research Fund for the Control of Infectious Diseases of the Health, Welfare and Food Bureau; the Providence Foundation Limited, in memory of the late Dr. Lui Hac Minh; and the Consultancy Service for Enhancing Laboratory Surveillance of Emerging Infectious Disease for the HKSAR Department of Health. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Introduction

More than 70% of the emerging infectious disease agents are caused by microbes jumping from animals into human. This has been well exemplified by the highly fatal human infection due to avian influenza A H5N1 in 1997 [1]. The outbreak of severe acute respiratory syndrome (SARS) caused by a novel coronavirus in 2003 [2], confirmed again that microbes can jump species from animals to humans with unpredictable consequence. The human SARS coronavirus was traced to caged civets in the market [3], and later Chinese horseshoe bat, Rhinolophus sinicus, was suggested to be a likely reservoir of SARS coronavirus [4]. Bats are ideal incubators for new emerging infectious agents as they are mammals which roosted together and can fly over vast geographical distance [5]. This has reignited the interest in seeking for new bat viruses including many bat coronaviruses and the recent discovery of bat influenza virus [6]. Besides the SARS coronavirus, viruses in bats often infect human through intermediate hosts such as horses for Hendra virus, pigs for Nipah virus, and chimpanzees for Ebola virus [5]. It is therefore important to catalogue as comprehensively as possible the animal viruses present in wild life especially the bats and birds, the food animals such as pigs and cattles, the pet animals such as cats and dogs, and monkeys which are phylogenetically close to humans. Using consensus primer polymerase chain reaction (PCR) screening, we have been able to discover relatively closely related species of virus in many different animals [4], [7][23]. However more distant or novel families of virus can only be found by metagnenomics using deep sequencing with the newer generation sequencers [24],[25]. We report in this paper the discovery and characterization of a novel bat papillomavirus (PV) from rectal swab samples randomly collected from asymptomatic wild, food and pet animals using a metagenomic approach.

Materials and MethodsSample collection

This study was performed in strict accordance with local ordinance and the recommendations by the Committee on the Use of Live Animals in Teaching and Research (CULATR) at the University of Hong Kong. The sampling of live animals were approved under permit no. 1048-05 (bats and monkeys) and 2284-10 (stray dogs and cats). All sampling were performed by licensed veterinarians, and anesthesia was given where appropriate; every effort was made to minimize suffering.

Sample collection was carried out in 2006–2007, and was approved by and performed in collaboration with the Department of Agriculture, Fisheries and Conservation (AFCD) of the Hong Kong Special Administrative Region (HKSAR). Collection of animal samples was performed by authorized staff members from AFCD of the HKSAR Government under the supervision of licensed veterinarian from AFCD, HKSAR (http://www.afcd.gov.hk/english/quarantine/qua_awc/qua_awc_leg/qua_awc_leg_dogs/qua_awc_leg_dogs.html).

10.1371/journal.pone.0043986.t001Size and position of predicted ORFs and NCR of MscPV1 and the predicted molecular masses of the translated proteins.
ORFPositionLengthMolecular mass (kDa)pI
ntaa
E61–42642614116.59.84
E7416–691276919.84.99
E1678–2687201066975.95.62
E22626–3768114338043.49.73
E43227–352930310011.56.42
L23946–5538159353057.56.36
L15546–7051150650156.66.84
NCR7052–7531480159NANA

A total of 96 rectal swabs were collected into viral transport medium from 10 types of bats including Rhinolophus sinicus (n = 10), Rhinolophus affinus (n = 10), Hipposideros pomona (n = 16), Miniopterus pusillus (n = 10), Miniopterus schreibersii (n = 10), Pipistrellus abramus (n = 10), Pipstrellus spp (n = 9), Myotis ricketti (n = 8), Myotis chinensis (n = 8), and Nyctalus noctula (n = 5). These bats were captured and sampled at 20 different locations in rural areas of the HKSAR, including water tunnels, abandoned mines, sea caves, and forested areas during a 1-year period. Bats were caught by nets during routine conservation procedures by AFCD, HKSAR. Collection of specimens was performed by an authorized veterinarian at the AFCD. Rectal swabs were collected from bats with medium-moistened cotton swab immediately immersed in viral transport medium. These bats were released after sample collection. The samples were collected for a routine surveillance study by AFCD.

The rectal swabs of 9 pigs and 9 cattles were collected and put into viral transport medium in a slaughter house in the New Territories, Hong Kong (Sheung Shui Slaughterhouse), a facility owned and operated by the HKSAR Government. The AFCD is a government department legitimately allowed to perform their duties in collaboration with other departments. The pigs and cattle had been previously slaughtered at the slaughterhouse. Samples from the carcasses were collected for a surveillance study by AFCD staff with departmental authorization.

Rectal swabs of 11 stray dogs, 11 stray cats and 20 monkeys were collected under anaesthesia. The least traumatic techniques were employed for the collection of samples. To minimize sufferings and injury, anaesthesia for restraining were carried out when necessary. The strays dogs and strays cats were kept at AFCD under standard facilities. The stray dogs and cats were euthanized after samplings, as part of the routine procedure of the AFCD. Wild monkeys were caught, temporarily kept in cages for less than one day, sampled and released, also as part of the routine procedure of the AFCD. A licensed AFCD veterinarian was responsible for assessing the well-being of animals for few hours and ensured that they are clinically normal before it was released back to nature. Procedures requiring institutional approval were approved by the Committee on the Use of Live Animals in Teaching and Research (CULATR) at the University of Hong Kong, permit numbers 1048-05 and 2284-10.

10.1371/journal.pone.0043986.t002MscPV1 nucleotide and amino acid identities with members of the genera <italic>Kappapapillomavirus, Lambdapapillomavirus, Mupapillomavirus, Nupapillomavirus, Sigmapapillomavirus,</italic> and <italic>Psipapillomavirus.</italic>
VirusaGenBank accession noGenusL1L2E1E2
ntaantaantaantaa
CPV6FJ492744Lambda60.058.645.735.454.944.448.040.3
CPV1D55633Lambda59.257.447.837.653.643.248.741.3
PlpPV1AY904724Lambda59.259.545.437.352.742.447.839.4
PlPV1AY763115Lambda59.059.146.436.954.546.151.045.3
EdPV1AY684126Sigma58.958.344.532.049.937.445.534.6
FdPV1AF480454Lambda58.959.145.437.553.343.448.338.9
HPV1V01116Mu58.856.545.435.453.443.548.339.2
UuPV1DQ180494Lambda58.858.445.338.252.542.947.140.5
SfPV1K02708Kappa58.755.943.731.851.540.844.735.6
LrPV1AY904722Lambda58.659.045.537.752.643.050.541.4
PcPV1AY904723Lambda58.658.445.036.453.546.047.839.2
HPV41X56147Nu58.156.443.834.148.336.444.132.3
HPV63X70828Mu57.657.946.336.052.342.150.140.7
OcPV1AF227240Kappa57.555.643.532.752.243.548.538.3
RaPV1DQ366842Psi52.444.840.027.648.738.345.331.7

CPV6, Canine papillomavirus 6; CPV1, Canine oral papillomavirus; PlpPV1, Panthera leo persica papillomavirus type 1; PlPV1, Procyon lotor papillomavirus type 1; EdPV1, Erethizon dorsatum papillomavirus type 1; FdPV1, Felis domesticus papillomavirus type 1; HPV1, Human papillomavirus type 1a; UuPV1, Uncia uncia papillomavirus type 1; SfPV1, Sylvilagus floridanus papillomavirus type 1; LrPV1, Lynx rufus papillomavirus type 1; PcPV1, Puma concolor papillomavirus type 1; HPV41, Human papillomavirus type 41; HPV63, Human papillomavirus type 63; OcPV1, Oryctolagus cuniculus papillomavirus type 1; RaPV1, Rousettus aegyptiacus papillomavirus type 1.

Sample preparation

The viral transport medium (in which the rectal swab specimens were immersed) were pooled, 100 μl each, and centrifuged at 10000×g for 5 min. The supernatant was then filtered through a 0.22 μm filter (Millipore). The filtrates were treated with DNaseI (Roche) and RNaseA (QIAGEN) to remove any extracellular nucleic acids that remained. Total RNA and DNA from the samples were extracted using the QIAamp Viral RNA Mini Kit (QIAGEN) and QIAamp DNA Mini Kit (QIAGEN), respectively. For the total RNA sample obtained, reverse transcription was performed using SuperScript III reverse transcriptase (Invitrogen) and random hexamers (Invitrogen) following the manufacturer's protocol. The cDNA and the previously extracted DNA were separately amplified using the Rapisome pWGA kit (Biohelix).

454 sequencing

The amplified DNA was used as a template for GS FLX analysis (Roche/454 Life Sciences) on one-quarter of a PicoTiterPlate according to manufacturer's instruction manual. The purified DNA level was determined by Nanodrop (Thermo Scientific). A total of 120 µg of DNA from the library was run on a 2% agarose gel, yielding a DNA smear. DNA ranging in size from 500 to 1,000 bp was cut from the gel and purified using the QIAquick Gel Extraction Kit (Qiagen). The extremities of the DNA fragments were then polished using T4 polynucleotide kinase. The Roche/454 adaptors were then ligated, and small DNA fragments removed and loaded on the machine according to the manufacturer's protocol (GS FLX Titanium General Library Preparation Kit, Roche).

Analysis of sequence reads

Sequences were trimmed based on quality score of 99.9% and any sequences less than 30 bp long were deleted. Duplicate sequences were discarded using 454 duplicate clustering workflows (sequence identity threshold of 0.96) based on the CDHIT program [26]. Remaining sequences were compared to a database containing 3,959 complete eukaryotic viral genomes (http://www.ncbi.nlm.nih.gov/genomes/VIRUSES/viruses.html) and the non-redundant protein sequences (nr) database from NCBI (http://www.ncbi.nlm.nih.gov) using tBLASTx and BLASTx, respectively, with an E-value cutoff of 10−5 [27]. BLAST results were parsed to save the best hits for each sequence. The best-hit sequences were individually annotated to note the sources of the matching sequences (eukaryotic virus, phage, bacteria and eukaryotes). Sequences were also analyzed using a metagenomic annotation tool, MEGAN version 4.50.6 to assign each sequence into different taxa present in the metagenomic sequences using the NCBI taxonomic database [28]. All unmapped reads were de novo assembled separately using MIRA to identify previously undetected virus [29].

<italic>De novo</italic> metagenomic assembly

De novo assembly of the metagenome was performed using MIRA to confirm isolation of viral genomes using an assembly option with minimum read length of 80 and base default quality of 10 [29]. There were 1,283 contigs ranging in size from 91 to 3,722 bp for human samples, and 1,960 contigs from 116 to 11,170 bp for animal samples. Contigs were compared to the database containing 3,959 complete eukaryotic viral genomes and the nr database from NCBI using tBLASTx and BLASTx, respectively, with an E-value cutoff of 10−5 to assign taxonomy.

Confirmation of the presence and completion of the assembled genome in the original specimen

After the incomplete genome of this novel virus was assembled from the metagenomic dataset, specific primers were designed to fill the gaps for the completion of this viral genome (primers available on request). For the confirmation of the host specimen containing this novel virus, DNA of each original specimen was subjected to PCR by using the primers specific for the L1 gene amplifying a 444 bp fragment between positions 483 to 926 (forward primer LPW11859 5′-GGCTCTCGGTGAGCACT-3′ and reverse primer LPW11861 5′-CAGTAAGGTCTGTTGAACAGTT-3′). The PCR mixture consisted of DNA template, PCR Buffer II at 1× (Applied Biosystems), 2 mM MgCl2, 200 μM of each dNTPs and 0.625 U AmpliTaq Gold DNA polymerase (Applied Biosystems). The mixtures were amplified in thermal cycler 9700 (Applied Biosystems), with a hot start of 95°C for 5 min, followed by 40 cycles of 95°C for 1 min, 55°C for 1 min and 72°C for 1 min and a final extension at 72°C for 10 min. PCR product was gel-purified using the QIAquick gel extraction kit (Qiagen). Both strands of the PCR products were sequenced with an ABI Prism 3700xl Genetic Analyser (Applied Biosystems) by using the PCR primers.

Distance measurements and phylogenetic analysis

The nucleotide global multiple sequence alignments were constructed for different open reading frames (ORFs) with 214 PVs based on the corresponding amino acid alignment using MUSCLE v3.7 [30] implemented in Seaview v4.1 as described previously [31], [32]. The pairwise identity values from nucleotides and proteins were calculated using MEGA5 [33]. Only the PV core early (E) ORFs E1 and E2 and the late (L) ORFs L1 and L2 were included as only these ORFs are ubiquitous present in all characterized PVs. L1 nucleotide sequences of MscPV1 and 78 PVs with complete genomes, representing all presently classified genera, were used for phylogenetic analysis. Maximum likelihood trees were constructed using PhyML with GTR+I+G model [34]. Modelgenerator was used to obtain the model for the likelihood analysis [35].

Genome analysis

Putative ORFs were predicted using ORF Finder and then searched for similarities with other proteins using BLASTP. Theoretical isoelectric points and molecular masses were estimated using Compute pI/Mw (http://web.expasy.org/compute_pi/). Proteins were analyzed for unique domains with InterProScan [36].

Prevalence of MscPV1 in bats

Prevalence of MscPV1 in M. schreibersii was further investigated by PCR screening of 419 additional samples (mouth swabs [n = 210], rectal swab [n = 127], anal swabs [n = 2], and urine samples [n = 80]) obtained from 257 bats using primers specific for the L1 gene amplifying a 444 bp fragment between positions 483 to 926 (forward primer LPW11859 5′-GGCTCTCGGTGAGCACT-3′ and reverse primer LPW11861 5′-CAGTAAGGTCTGTTGAACAGTT-3′).

Nucleotide sequence accession number

The nucleotide sequence of the genome of MscPV1 has been lodged within the GenBank sequence database under accession no. JQ692938.

ResultsIdentification of a novel papillomavirus

Approximately 3% of the sequence reads generated from these animal samples were assigned to eukaryotic viral sequences by the BLAST nr protein database. The majority of the viral-like sequences were similar to single-stranded, negative-sense, circular DNA viruses, with the largest proportion of the sequences showing homology to porcine circovirus. The next large group of the sequences matched to another member of the Circoviridae family, torque teno virus, including torque teno felis virus, torque teno sus virus 1 and torque teno canis virus. The remaining viral-like sequences shared homology to canary circovirus, anellovirus and densovirus in which densovirus is a linear single-stranded DNA virus. In addition, many sequences were categorized as phage-related genes (Fig. 1). The majority of these sequences were related to porcine circovirus in animal samples. Twenty two sequence reads and one contig in animal sample were related to PVs with amino acid identity ranged from 42% to 73%. These hits cover about 70% of the viral genome, which we named Miniopterus schreibersii papillomavirus type 1 (MscPV1), since the sample was isolated from a Common Bent-wing Bat (M. schreibersii). This bat is a female adult bat collected on 29 December 2006 in Tung Tsz, Hong Kong. By connecting gaps between sequenced viral fragments based on PV sequences, the complete genome of the novel PV was acquired.

10.1371/journal.pone.0043986.g001MEGAN tree with taxonomic assignments.

The distribution of the sequence reads through blastx analysis against the nr database. Size of circles located next to taxa are proportional to the total number of reads identified. Not assigned contains those reads that are not assigned by the least common ancestor algorithm. No hits contains those reads that did not return any significant alignments to the nr database.

Characterization of MscPV1 complete genome

The complete genome of MscPV1 was 7,531 bp in length with a G+C content of 45.8%. The MscPV1 genome contains the typical PV ORFs, coding for five putative early proteins (E6, E7, E1, E2, E4), and two putative late capsid proteins (L2 and L1) (Fig. 2 and Table 1).

10.1371/journal.pone.0043986.g002Circular and linear genome maps of <italic>Miniopterus schreibersii</italic> papillomavirus type 1 (MscPV1).

Characteristic features of the long control region of MscPV1, showing genomic locations of E2 binding sites (bold), polyadenylation sites (underlined), and TATA box (boxed).

The MscPV1 E6 contains two conserved zinc binding domains (CXXCX29CXXC), separated by 36 amino acids, whereas the MscPV1 E7 contained one slightly modified domain (CXXCX30CXXC), but no retinoblastoma tumour suppressor (pRB)-binding domain (LXCXE) [37]. The E1 ORF codes for the largest MscPV1 protein (669 aa), and contains the conserved ATP-binding site of the ATP-dependent helicase (GXXXXGK(T/S)) [38]. This sequence is GPPDTGKS in MscPV1. The E2 protein has the typical C-terminal DNA-binding domain and the N-terminal transactivation domain [39], [40]. The MscPV1 E4 gene is located in the E region and overlaps with E2 but is transcribed in a different reading frame. An LLXLL motif is found at the N-terminus of viral E4 [41]. Downstream from the leucine-rich region is a proline-rich region. PV E4 proteins usually have high proline content (15–20% on average), MscPV1 E4 protein also has the typical high proline content (15 proline residues out of 106 aa).

Both L1 and L2 contain a series of arginine and lysine residues at their carboxy termini, likely to function as a nuclear localization signal. The long control region (NCR) usually contains several regulators of the PV replication. In MscPV1, the NCR is 481 bp and demonstrates an E1-binding site (TGATTGTTGTAAACTAC) flanked by two typical palindromic E2-binding sites (ACCN6GGT) [42]. At its 5′ end, the NCR also contains one polyadenylation site (AATAAA) which is necessary for the processing of the L1 and L2 capsid mRNA transcript [43]. In the 3′ end, the MscPV1 NCR contains a classical TATA box (TATAAA) of the E6 promotor, located 26 nucleotides upstream of the E6 start codon (Fig. 2).

Phylogenetic analysis and sequence similarity to other papillomaviruses

Phylogenetic analysis confirmed that MscPV1 forms a genetic lineage that is distinct from the previously reported PVs with complete genome (Fig. 3). Comparison of L1 gene showed that MscPV1 had 60% nucleotide and 58.6% amino acid identity to the closest related PV, Canine papillomavirus 6 (Table 2). MscPV1 also shared only 52.4% nucleotide identity to another PV isolated from an Egyptian fruit bat (Rousettus aegyptiacus) (Table 2) [44]. MscPV1 cannot be placed in one of the existing genera, it therefore represents the first member of a novel PV genus, Dyolambda-papillomavirus, according to the classification criteria [31], [45].

10.1371/journal.pone.0043986.g003Maximum likelihood phylogenetic tree of the L1 nucleotide sequences of 79 PVs.

The PV genus of each strain is indicated. PVs with putative PV genera that are currently unclassified are marked by asterisks. The PV discovered in this study is shown in bold. Scale bar indicates 0.2 inferred substitutions per site. AaPV, Alces alces papillomavirus; BpPV, Bettongia penicillata papillomavirus; BPV, Bovine papillomavirus; CcaPV, Capreolus capreolus papillomavirus; CcPV, Caretta caretta papillomavirus; CgPV, Colobus guereza papillomavirus; ChPV, Capra hircus papillomavirus; CPV, Canine papillomavirus; EcPV, Equus caballus papillomavirus; EdPV, Erethizon dorsatum papillomavirus; EePV, Erinaceus europaeus papillomavirus; FcPV, Fringilla coelebs papillomavirus; FdPV, Felis domesticus papillomavirus; FlPV, Francolinus leucoscepus papillomavirus; HPV, Human papillomavirus; LrPV, Lynx rufus papillomavirus; MaPV, Mesocricetus auratus papillomavirus; MfPV, Macaca fascicularis papillomavirus; MmiPV, Micromys minutus papillomavirus; MmPV, Macaca mulatta papillomavirus; MnPV, Mastomys natalensis papillomavirus; MscPV, Miniopterus schreibersii papillomavirus; MsPV, Morelia spilota spilota papillomavirus; OaPV, Ovis aries papillomavirus; OcPV, Oryctolagus cuniculus papillomavirus; OvPV, Odocoileus virginianus papillomavirus; PcPV, Puma concolor papillomavirus; PePV, Psittacus erithacus timneh papillomavirus; PlpPV, Panthera leo persica papillomavirus; PlPV, Procyon lotor papillomavirus; PpPV, Pygmy chimpanzee papillomavirus; PsPV, Phocoena spinipinnis papillomavirus; RaPV, Rousettus aegyptiacus papillomavirus; SfPV, Sylvilagus floridanus papillomavirus; SsPV, Sus scrofa papillomavirus; TmPV, Trichechus manatus latirostris papillomavirus; TtPV, Tursiops truncatus papillomavirus; UmPV, Ursus maritimus papillomavirus; UuPV, Uncia uncia papillomavirus; ZcPV, Zalophus californianus papillomavirus.

Prevalence of MscPV1 in bats

None of the 419 samples from 257 M. Schreibersii bats screened by PCR was positive.

Discussion

Virus discovery has traditionally been done by phenotypic techniques such as animal inoculation or chick embryo inoculation which are later replaced by tissue culture. With major advance in molecular and sequencing technology, many viruses that may not adapt to grow in tissue cultures were discovered by PCR and sequencing in various formats such as consensus primer PCR with or without hybridization on microarray, rolling circle amplification for virus with circular genome, representational difference analysis by subtractive hybridization, sequence independent PCR amplification with shotgun sequencing. The advent of high-throughput sequencing has allowed the discovery of many novel animal viruses such as novel species of porcine circoviruses, astroviruses and bocaviruses [46], [47], novel sapoviruses, noroviruses, dependoviruses in sealions [48], novel kobuvirus and sapovirus in diarrheal dogs [49], novel hepacivirus in dogs affected by outbreak of respiratory illness [50], novel anellovirus in sea seals [51], novel astrovirus in brain tissue of mink suffering from shaking mink syndrome [52] and many other virus families in human [53], turkey [54], bat guano [55], rodent excreta [56] and insects [57], [58]. Metagenomics has also led to the study of the viral diversity and community in hosts and the association of virus and disease [59][61].

In this study, we report the second bat PV, MscPV1, with only 52.4% nucleotide identity to the Egyptian fruit bat (Rousettus aegyptiacus) papillomavirus RaPV1. Compared to other PVs, the highest nucleotide and amino acid identities, from CPV6, were only 60% and 58.6%. According to the published classification criteria, MscPV1 should be designated the first member of a novel PV genus, Dyolambda-papillomavirus [31], [45].

Comparing MscPV1 with the phylogenetically closely-related PVs, namely, RaPV1, HPV41, EdPV1 and CPV6, all of MscPV1, RaPV1, EdPV1 and CPV4 contain the typical PV ORFs, coding for five putative early proteins (E6, E7, E1, E2, E4), and two putative late capsid proteins (L2 and L1). The genome of HPV41 consists of an additional E5 ORF located between E4 and L2 ORFs and three additional short ORFs, X, Y and Z downstream of L1 [62]. The E5 ORF, which is absent in MscPV1, exists in genital HPVs and in the BPV-1 related fibropaillomaviruses, codes for the E5 protein, which is associated with transformation of host cells and carcinogenesis [62], [63]. The predicted E7 protein of MscPV1 contains a modified zinc-binding domain with 30 amino acids (X30) between the two instances of CXXC. This nonclassical motif was also identified in HPV41 and RaPV1 as well as BPV6, CPV2, CPV7, CcaPV1, HPV4, HPV65, HPV95 and HPV116. The E7 of EdPV1 exhibits the classical CXXCX29CXXC motif, whereas CPV6 has the X28 modified motif. The E1 of MscPV1 contains the conserved ATP-binding site of the ATP-dependent helicase (GPPDTGKS), which is identical to that of RaPV1; in comparison, the motif is GPPNTGKS in EdPV1 and CPV6, and GPSDTGKS in HPV41. This sequence conservation is not unexpected, given the drastic decrease of ATPase activity upon mutation of just the first proline or the lysine residue of the motif demonstrated by a mechanistic study [64]. In the NCR of both MscPV1 and RaPV1, two copies of the 12-basepair E2 protein-binding motif ACCN6GGT are found. The genomes of HPV41 and EdPV1, notably, do not contain this consensus nucleotide sequence; their E2 binding sites are represented by the sequences ACCN6GTT, AACN6GGT, and AACN6GTT [62], [65].

Family Papillomaviridae is a large family of small, non-enveloped, double stranded DNA viruses which infect cutaneous and mucosal epithelium. Given the association between PV and cancers in humans and other animals, it is not surprising that the first bat PV was found in the basosquamous carcinoma of a fruit bat. PVs are stable and slow-evolving viruses, with an estimated mutation rate of 0.73 to 1.2×10−8 nucleotide substitutions per base per year [66], [67]. No genomic recombination has ever been documented. Novel PVs, therefore, have been believed to descend from the slow accumulation of point mutations and different ancient PV lineages have possibly co-evolved and co-speciated with their vertebrate host species [68], [69]. Nonetheless, a more recent study showed that recombination of PV contributes significantly to the evolution of PV, and the role of host transfer cannot be neglected. Notably, both were observed in the HPV41-EdPV1 clade [70]. This is concordant with the current study, in which nucleotide and amino acid sequence analysis demonstrated a higher degree of similarity between our novel MscPV1 and EdPV1 in North American porcupine, instead of another bat PV, RaPV1. We note, however, our findings may only represent an additional exception and do not refute the generalization that PVs evolve mainly by co-evolution with hosts instead of interspecies transmission, recombination or other horizontal genetic transfer events. While current evidence suggests that host shift may have contributed to the emergence of this lineage of bat PV, as evidenced by the different clades in which the MscPV1 and we reckon that there may be bat PV lineages that are yet to be discovered, and such “missing-links” may enable the evolutionary history of MscPV1 to be put into context, i.e. whether the MscPV1 has undergone unexpected, rapid divergent evolution or indeed represents a lineage of bat PVs that has arisen from host shift. Although divergence time estimation using the L1 nucleotide sequence of MscPV1 demonstrated a later divergence compared to the host species (Fig. S1), it differed from the estimate based on E1 nucleotide sequences (Fig. S2) and does not in itself substantiate the host transfer (Supporting Information S1). In the current study, the more convincing evidence seemingly still resides in the sister taxa status of MscPV1 and HPV41/EdPV1. The relative importance of virus-host co-divergence and interspecies transmission in driving the genomic evolution of PV remains to be debated [71].

Unlike the first report of bat PV, no obvious skin or mucosal lesion was noted by the attending veterinarian and the bat was probably latently infected with PV. M. schreibersii is a cave-dwelling bat with body weight ranged from 11 to 18 g. It roosts in abandoned mines, water tunnels, drainage and weep holes of the water catchment. According to the baseline surveys of the AFCD of Hong Kong, M. schreibersii is considered as common and widespread throughout Hong Kong countryside with a colony size from 50 individuals to several hundreds, often associated with M. pusillus, Myotis pilosus and M. chinensis in their roosting sites in summer. In overseas studies, it is reported as a migratory bat species which may travel a fairly long distance in spring to find their breeding sites [72]. Given the possibilities of asymptomatic carriage and long-distance, interspecies transmission, further studies are warranted to elucidate the evolutionary origin and epidemiology of this newly proposed genus of bat PV.

Supporting Information

Estimation of the time to the most recent common ancestor for MscPV1 using L1. The maximum likelihood tree constructed by PhyML using L1 were used to estimate the divergence times in MEGA5. Virus name abbreviations are the same as those in the Fig. 3 legend. MscPV1 was bolded.

(TIF)

Click here for additional data file.

Estimation of the time to the most recent common ancestor for MscPV1 using E1. The maximum likelihood tree constructed by PhyML using E1 were used to estimate the divergence times in MEGA5. Virus name abbreviations are the same as those in the Fig. 3 legend. MscPV1 was bolded.

(TIF)

Click here for additional data file.

(DOC)

Click here for additional data file.

We thank Alan Chi-Kong Wong, Siu-Fai Leung, Chik-Chuen Lay, Thomas Sit, K. F. Chan, Chung-Tong Shek, Cynthia S. M. Chan, Joseph W. K. So, Michelle L. Yeung, Byung Mo Hwang, Suet Yee Ng, Patrick I. T. Lau, and Steven D. Benton from HKSAR Department of Agriculture, Fisheries, and Conservation (AFCD) for facilitation and support and members of Animal Management Centres of AFCD. We are grateful for the generous support of Carol Yu, Richard Yu, Hui Hoy, and Hui Ming with the genomic sequencing platform.

ReferencesYuenKY, ChanPK, PeirisM, TsangDN, QueTL, et al (1998) Clinical features and rapid viral diagnosis of human disease associated with avian influenza A H5N1 virus. Lancet351: 4674719482437PeirisJS, LaiST, PoonLL, GuanY, YamLY, et al (2003) Coronavirus as a possible cause of severe acute respiratory syndrome. Lancet361: 1319132512711465GuanY, ZhengBJ, HeYQ, LiuXL, ZhuangZX, et al (2003) Isolation and characterization of viruses related to the SARS coronavirus from animals in southern China. Science302: 27627812958366LauSK, WooPC, LiKS, HuangY, TsoiHW, et al (2005) Severe acute respiratory syndrome coronavirus-like virus in Chinese horseshoe bats. Proc Natl Acad Sci U S A102: 140401404516169905WongS, LauS, WooP, YuenKY (2007) Bats as a continuing source of emerging infections in humans. Rev Med Virol17: 679117042030TongS, LiY, RivaillerP, ConrardyC, CastilloDA, et al (2012) A distinct lineage of influenza A virus from bats. Proc Natl Acad Sci U S A109: 4269427422371588Woo PC, Lau SK, Lam CS, Lau CC, Tsang AK, et al (2012) Discovery of seven novel mammalian and avian coronaviruses in Deltacoronavirus supports bat coronaviruses as the gene source of Alphacoronavirus and Betacoronavirus and avian coronaviruses as the gene source of Gammacoronavirus and Deltacoronavirus J Virol.WooPC, LauSK, ChuCM, ChanKH, TsoiHW, et al (2005) Characterization and complete genome sequence of a novel coronavirus, coronavirus HKU1, from patients with pneumonia. J Virol79: 88489515613317LauSK, WooPC, LiKS, HuangY, WangM, et al (2007) Complete genome sequence of bat coronavirus HKU2 from Chinese horseshoe bats revealed a much smaller spike gene with a different evolutionary lineage from the rest of the genome. Virology367: 42843917617433LauSK, PoonRW, WongBH, WangM, HuangY, et al (2010) Coexistence of different genotypes in the same bat and serological characterization of Rousettus bat coronavirus HKU9 belonging to a novel Betacoronavirus subgroup. J Virol84: 113851139420702646LauSK, LiKS, HuangY, ShekCT, TseH, et al (2010) Ecoepidemiology and complete genome comparison of different strains of severe acute respiratory syndrome-related Rhinolophus bat coronavirus in China reveal bats as a reservoir for acute, self-limiting infection that allows recombination events. J Virol84: 2808281920071579WooPC, LauSK, LiKS, PoonRW, WongBH, et al (2006) Molecular diversity of coronaviruses in bats. Virology351: 18018716647731WooPC, WangM, LauSK, XuH, PoonRW, et al (2007) Comparative analysis of twelve genomes of three novel group 2c and group 2d coronaviruses reveals unique group and subgroup features. J Virol81: 1574158517121802LauSK, WooPC, LaiKK, HuangY, YipCC, et al (2011) Complete genome analysis of three novel picornaviruses from diverse bat species. J Virol85: 8819882821697464LauSK, WooPC, WongBH, WongAY, TsoiHW, et al (2010) Identification and complete genome analysis of three novel paramyxoviruses, Tuhoko virus 1, 2 and 3, in fruit bats from China. Virology404: 10611620537670LauSK, WooPC, TseH, FuCT, AuWK, et al (2008) Identification of novel porcine and bovine parvoviruses closely related to human parvovirus 4. J Gen Virol89: 1840184818632954TseH, TsoiHW, TengJL, ChenXC, LiuH, et al (2011) Discovery and genomic characterization of a novel ovine partetravirus and a new genotype of bovine partetravirus. PLoS One6: e2561921980506WooPC, LauSK, LamCS, LaiKK, HuangY, et al (2009) Comparative analysis of complete genome sequences of three avian coronaviruses reveals a novel group 3c coronavirus. J Virol83: 90891718971277LauSK, WooPC, YipCC, ChoiGK, WuY, et al (2012) Identification of a novel feline picornavirus from the domestic cat. J Virol86: 39540522031936WooPC, LauSK, HuangY, LamCS, PoonRW, et al (2010) Comparative analysis of six genome sequences of three novel picornaviruses, turdiviruses 1, 2 and 3, in dead wild birds, and proposal of two novel genera, Orthoturdivirus and Paraturdivirus, in the family Picornaviridae. J Gen Virol91: 2433244820554801WooPC, LauSK, WongBH, WongAY, PoonRW, et al (2011) Complete genome sequence of a novel paramyxovirus, Tailam virus, discovered in Sikkim rats. J Virol85: 134731347422106385TseH, ChanWM, TsoiHW, FanRY, LauCC, et al (2011) Rediscovery and genomic characterization of bovine astroviruses. J Gen Virol92: 1888189821508185LauSK, WooPC, YipCC, LiKS, FuCT, et al (2011) Co-existence of multiple strains of two novel porcine bocaviruses in the same pig, a previously undescribed phenomenon in members of the family Parvoviridae, and evidence for inter- and intra-host genetic diversity and recombination. J Gen Virol92: 2047205921632566EndyTP, RochfordR, YuenKY, LeiHY (2011) Emerging infectious diseases as a global health threat. Introduction. Exp Biol Med (Maywood)236: 89789821821661FinkbeinerSR, AllredAF, TarrPI, KleinEJ, KirkwoodCD, et al (2008) Metagenomic analysis of human diarrhea: viral detection and discovery. PLoS Pathog4: e100001118398449LiW, GodzikA (2006) Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics22: 1658165916731699AltschulSF, GishW, MillerW, MyersEW, LipmanDJ (1990) Basic local alignment search tool. J Mol Biol215: 4034102231712HusonDH, AuchAF, QiJ, SchusterSC (2007) MEGAN analysis of metagenomic data. Genome Res17: 37738617255551ChevreuxB, WetterT, SuhaiS (1999) Genome Sequence Assembly Using Trace Signals and Additional Sequence Information. Computer Science and Biology: Proceedings of the German Conference on Bioinformatics (GCB)99: 4556EdgarRC (2004) MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res32: 1792179715034147BernardHU, BurkRD, ChenZ, van DoorslaerK, HausenH, et al (2010) Classification of papillomaviruses (PVs) based on 189 PV types and proposal of taxonomic amendments. Virology401: 707920206957GouyM, GuindonS, GascuelO (2010) SeaView version 4: A multiplatform graphical user interface for sequence alignment and phylogenetic tree building. Mol Biol Evol27: 22122419854763TamuraK, PetersonD, PetersonN, StecherG, NeiM, et al (2011) MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol28: 2731273921546353GuindonS, DufayardJF, LefortV, AnisimovaM, HordijkW, et al (2010) New algorithms and methods to estimate maximum–likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol59: 30732120525638KeaneTM, CreeveyCJ, PentonyMM, NaughtonTJ, McLnerneyJO (2006) Assessment of methods for amino acid matrix selection and their use on empirical data shows that ad hoc assumptions for choice of matrix are not justified. BMC Evol Biol6: 2916563161ZdobnovEM, ApweilerR (2001) InterProScan-an integration platform for the signature-recognition methods in InterPro. Bioinformatics17: 84784811590104DahiyaA, GavinMR, LuoRX, DeanDC (2000) Role of the LXCXE binding site in Rb function. Mol Cell Biol20: 6799680510958676TitoloS, PelletierA, SauveF, BraultK, WardropE, et al (1999) Role of the ATP-binding domain of the human papillomavirus type 11 E1 helicase in E2-dependent binding to the origin. J Virol73: 5282529310364274HegdeRS, GrossmanSR, LaiminsLA, SiglerPB (1992) Crystal structure at 1.7 A of the bovine papillomavirus-1 E2 DNA-binding domain bound to its DNA target. Nature359: 5055121328886HegdeRS (2002) The papillomavirus E2 proteins: structure, function, and biology. Annu Rev Biophys Biomol Struct31: 34336011988474RobertsS, AshmoleI, GibsonLJ, RookesSM, BartonGJ, et al (1994) Mutational analysis of human papillomavirus E4 proteins: identification of structural features important in the formation of cytoplasmic E4/cytokeratin networks in epithelial cells. J Virol68: 643264457521917LiR, KnightJ, BreamG, StenlundA, BotchanM (1989) Specific recognition nucleotides and their DNA context determine the affinity of E2 protein for 17 binding sites in the BPV-1 genome. Genes Dev3: 5105262542129BirnstielML, BusslingerM, StrubK (1985) Transcription termination and 3′ processing: the end is in site! Cell. 41: 349359RectorA, MostmansS, Van DoorslaerK, McKnightCA, MaesRK, et al (2006) Genetic characterization of the first chiropteran papillomavirus, isolated from a basosquamous carcinoma in an Egyptian fruit bat: the Rousettus aegyptiacus papillomavirus type 1. Vet Microbiol117: 26727516854536de VilliersEM, FauquetC, BrokerTR, BernardHU, zur HausenH (2004) Classification of papillomaviruses. Virology324: 172715183049BaylisSA, FinsterbuschT, BannertN, BlumelJ, MankertzA (2011) Analysis of porcine circovirus type 1 detected in Rotarix vaccine. Vaccine29: 69069721093497ShanT, LiL, SimmondsP, WangC, MoeserA, et al (2011) The fecal virome of pigs on a high-density farm. J Virol85: 116971170821900163LiL, ShanT, WangC, CoteC, KolmanJ, et al (2011) The fecal viral flora of California sea lions. J Virol85: 9909991721795334LiL, PesaventoPA, ShanT, LeuteneggerCM, WangC, et al (2011) Viruses in diarrhoeic dogs include novel kobuviruses and sapoviruses. J Gen Virol92: 2534254121775584KapoorA, SimmondsP, GeroldG, QaisarN, JainK, et al (2011) Characterization of a canine homolog of hepatitis C virus. Proc Natl Acad Sci U S A108: 116081161321610165NgTF, WheelerE, GreigD, WaltzekTB, GullandF, et al (2011) Metagenomic identification of a novel anellovirus in Pacific harbor seal (Phoca vitulina richardsii) lung samples and its detection in samples from multiple years. J Gen Virol92: 1318132321402596BlomstromAL, WidenF, HammerAS, BelakS, BergM (2010) Detection of a novel astrovirus in brain tissue of mink suffering from shaking mink syndrome by use of viral metagenomics. J Clin Microbiol48: 4392439620926705GreningerAL, ChenEC, SittlerT, ScheinermanA, RoubinianN, et al (2010) A metagenomic analysis of pandemic influenza A (2009 H1N1) infection in patients from North America. PLoS One5: e1338120976137DayJM, BallardLL, DukeMV, SchefflerBE, ZsakL (2010) Metagenomic analysis of the turkey gut RNA virus community. Virol J7: 31321073719LiL, VictoriaJG, WangC, JonesM, FellersGM, et al (2010) Bat guano virome: predominance of dietary viruses from insects and plants plus novel mammalian viruses. J Virol84: 6955696520463061PhanTG, KapusinszkyB, WangC, RoseRK, LiptonHL, et al (2011) The fecal viral flora of wild rodents. PLoS Pathog7: e100221821909269NgTF, WillnerDL, LimYW, SchmiederR, ChauB, et al (2011) Broad surveys of DNA viral diversity obtained through viral metagenomics of mosquitoes. PLoS One6: e2057921674005RosarioK, MarinovM, StaintonD, KrabergerS, WiltshireEJ, et al (2011) Dragonfly cyclovirus, a novel single-stranded DNA virus discovered in dragonflies (Odonata: Anisoptera). J Gen Virol92: 1302130821367985DonaldsonEF, HaskewAN, GatesJE, HuynhJ, MooreCJ, et al (2010) Metagenomic analysis of the viromes of three North American bat species: viral diversity among different bat species that share a common habitat. J Virol84: 130041301820926577SullivanPF, AllanderT, LysholmF, GohS, PerssonB, et al (2011) An unbiased metagenomic search for infectious agents using monozygotic twins discordant for chronic fatigue. BMC Microbiol11: 221194495TokarzR, FirthC, StreetC, Cox-FosterDL, LipkinWI (2011) Lack of evidence for an association between Iridovirus and colony collapse disorder. PLoS One6: e2184421738798HirtL, Hirsch-BehnamA, de VilliersEM (1991) Nucleotide sequence of human papillomavirus (HPV) type 41: an unusual HPV type without a typical E2 binding site consensus sequence. Virus Res18: 1791891645904DiMaioD, MattoonD (2001) Mechanisms of cell transformation by papillomavirus E5 proteins. Oncogene20: 7866787311753669WhitePW, PelletierA, BraultK, TitoloS, WelchnerE, et al (2001) Characterization of recombinant HPV6 and 11 E1 helicases: effect of ATP on the interaction of E1 with E2 and mapping of a minimal helicase domain. J Biol Chem276: 224262243811304544RectorA, TachezyR, Van DoorslaerK, MacNamaraT, BurkRD, et al (2005) Isolation and cloning of a papillomavirus from a North American porcupine by using multiply primed rolling-circle amplification: the Erethizon dorsatum papillomavirus type 1. Virology331: 44945615629787Van Ranst M, Sundberg JP, Burk RD (1995) Molecular Basis of Virus Evolution. Cambridge University Press: 455–476.TachezyR, DusonG, RectorA, JensonAB, SundbergJP, et al (2002) Cloning and genomic characterization of Felis domesticus papillomavirus type 1. Virology301: 31332112359433Van RanstM, FuseA, FitenP, BeukenE, PfisterH, et al (1992) Human papillomavirus type 13 and pygmy chimpanzee papillomavirus type 1: comparison of the genome organizations. Virology190: 5875961325697TachezyR, RectorA, HavelkovaM, WollantsE, FitenP, et al (2002) Avian papillomaviruses: the parrot Psittacus erithacus papillomavirus (PePV) genome has a unique organization of the early protein region and is phylogenetically related to the chaffinch papillomavirus. BMC Microbiol2: 1912110158ShahSD, DoorbarJ, GoldsteinRA (2010) Analysis of host-parasite incongruence in papillomavirus evolution using importance sampling. Mol Biol Evol27: 1301131420093429GottschlingM, GokerM, StamatakisA, Bininda-EmondsOR, NindlI, et al (2011) Quantifying the phylodynamic forces driving papillomavirus evolution. Mol Biol Evol28: 2101211321285031Kunz TH, Fenton MB (2003) Bat Ecology: University of Chicago Press. 779 p.