Identification and Characterization of Eleven Novel Human Gamma-papillomavirus Isolates from Healthy Skin, Found at Low Frequency in a Normal Population

Eleven novel human papillomavirus (HPV) types were isolated and characterized from healthy individuals in China. HPV163 belongs to the γ-1 species, HPV 164 and HPV 168 fit in the γ-8 species, HPV 165 and KC5 belongs to the γ-12 species, HPV 168 is closely allied with the γ-4 species, HPV 169 is closely related to the γ-11 species, and HPV 170 is related to the γ-12 species. In addition, HPV 161, HPV 162, and HPV 166 may form a new HPV species of the γ-PV genus. The prevalence of these HPV types in the normal population is low.


Introduction
Papillomavirus (PVs) are a distinct viral family, which infect the epithelia of vertebrates. Infections may be prolonged without clinical symptoms; however, they may induce neoplasia in certain sites [1,2]. PV infection is species restricted and is rarely transmitted across species. Human papillomavirus is one of the most important groups of the PV family, and it can be divided roughly into two tropism groups, cutaneotropic and mucosotropic HPV types, which infect skin and mucosal tissues respectively. Because lack of effective in vivo PV culture technology and PVs do not naturally induce robust antibody production in host, traditional virus classification is not appropriate for PVs. A typical novel PV type is identified and defined by cloning full-length PV genome whose L1 gene nucleotides sequences is at least 10% dissimilar from any other known PV type. The classification of PV types is based on construction of a phylogenetic tree using the L1 genes [3]. Within the papillomaviridae, 37 separate genera have been established, and five of them infect human beings. Five human papillomavirus genera including Alpha, Beta, Gamma, Mu, and Nu contain at least 150 HPV types that have been fully characterized. However, it is believed that there is potentially much greater number of unknown HPV types to be identified, and most of the undiscovered HPVs are expected to be cutaneotropic types.
Compared with mucosotropic HPV types, cutaneotropic HPV types exhibited very different characteristics not only in their nucleotide sequences but also in their biologic behavior. It is well known that approximately 13 high-risk mucosotropic HPV types can cause cervical cancer (IARC monographs), but none of the cutaneotropic HPV types has been confirmed to be associated with any malignant skin disease [4], which indicates that the cutaneotropic HPVs may exhibit lower pathogenicity. The prevalence and diversity of the cutaneotropic HPV group are both much higher than mucosotropic HPVs in a healthy population [5][6][7]. In addition, normal human skin harbors an array of papillomaviruses, and most of them are unknown previously [8]. Therefore there is great potential for discovery of novel HPV types from the normal skin of healthy individuals. Isolation and characterization novel HPV genomes will reveal the existence of distinct human papillomavirus (HPV) types and will provide basic information for understanding the biologic significance and the basis of HPV evolution. Many studies focus on identifying novel HPV from common warts, keratotic lesions or cutaneous squamous cell carcinoma [9][10][11][12][13][14], but few studies have investigated novel HPV types originated from healthy individuals.
In this study, full-length genomes of eleven novel HPV types were isolated, fully cloned, sequenced and characterized from healthy individuals in rural Anyang, Henan province, China. All of these eleven novel HPVs belong to the Gamma-papillomavirus genus, and show very low prevalence in this population.

Skin samples
Exfoliated skin cell samples used in this study were collected from individuals who participated in a perspective cohort study of esophageal cancer launched in rural An Yang, China [15]. For each individual, exfoliated cells were collected by scraping the palms of the hands using glass slides and then collecting the resultants material into 1 ml of 0.9% sodium chloride. Cells were centrifuged at 10,000 rpm for 10 min. Precipitates were collected and stored in -80°C before use. DNA was extracted using the spin-column method (Tiagen, a subsidiary of Qiagen in China), and was stored at -20°C before analyzed.

Isolation and characterization of novel HPV types
Novel HPV clones were initially identified by PCR amplification with FAP6085/64 primer pairs [16,17] producing a 377 bp L1 fragment. A total of 149 putative novel HPV types were identified in detecting skin HPV infection from a healthy population [18], and 11 of them were successfully whole genome sequenced and cloned into vectors.
Rolling circle amplification (RCA) was performed to amplify the circle virus genome DNA by using illustra TempliPhi 100 Amplification Kit (GE Healthcare, Amersham, UK), as previous described [19]. RCA products were diluted 20-100 fold as a PCR template. Type-specific primers were designed for each putative novel HPV types based on the FA amplified fragment sequences. Reverse PCR was performed using the NEB LongAmp Taq system (New England Biolabs). The 30 µl PCR reaction buffer consisted of 0.4 µM of each primer, 200 µM of each dNTP, 1x LongAmp buffer, 2.5 mM Mg 2 SO 4 and 2.5 U LongAmp Taq enzyme, and 1 µl of diluted RCA product template. PCR was initiated by a preheating step at 94°C for 1 min, followed by 40 cycles consisting of 94 °C for 20 s, 60°C for 40 s, and 65 °C for 7 min. The final elongation step was at 65°C for 10 min. The PCR amplicons were excised from an agarose gel and cloned into a pCR 2.1 TOPO TA vector (Invitrogen), according to the manufacturers' instructions, and sequenced with M13 primers.
To construct HPV clones containing whole genome sized fragment, other sets of type-specific primers were designed by using Primer3 software [20]. With this design, each HPV gene remained intact except in the LCR region ( Figure 1). For HPV 161-166, HPV 169, HPV 170, and KC5, type-specific longtemplate PCR was performed using the NEB LongAmp Taq system (New England Biolabs) with the same reaction system and conditions as described above, with primers as listed in Table S1. For Cloning HPV 167 and HPV 168, two overlap fragments were amplified by two pairs of primers (Table S2). PCR amplifying overlap fragments were performed in a 30 µl reaction buffer with 0.5µM primers, 200µM each of dNTP, 1x Phusion HF buffer, 0.5 U Phusion DNA polymerase (New England Biolabs), and 1 µl of diluted RCA product as a template. The conditions for high-fidelity PCR were: 98°C for 30s, and then 35 cycles at 98°C for 10s, 60°C for 15s, and 72°C for 4 min. The final extension was performed at 72°C for 10 min. Blunted PCR products were cloned by using the CloneJET TM PCR Cloning Kit (Fermentas) according to the manufacturer's instructions. Full sequencing by primer walking was carried out for all plasmid clones containing whole HPV genomes. Plasmids clones and sequences were submitted to the Reference Centre for Papillomaviruses, Heidelberg, Germany for official consideration. The complete genome sequences were submitted to the NCBI/GenBank database.

Protein Prediction
ORFs were predicted by using the ORFinder (http:// www.ncbi.nlm.nih.gov/gorf/gorf.ht-ml). E4 genes were identified by searching the proline-rich region (private communication with Dr. de Villiers, Dr. Eric Schulz, and Dr. Boštjan Kocjan) and confirmed by manual alignment of nucleotide and amino acid sequences to the homologous regions of related HPV types. Further genome structure analysis was predicted by directly searching the expected conserved protein or nucleotide sequences [21] with assistance of the CLC sequence view 6 program (http://www.clcbio.com/products/clc-sequenceviewer/).

Sequence alignment and Phylogenetic analysis
Multiple alignment of full-length L1 ORF of 11 novel HPV and 57 previously characterized HPV types obtained from the PaVE database (http://pave.niaid.nih.gov/#home) was performed by using the MEGA5 software package [22]. Fifty-seven previously knew HPVs include 35 published sequences of the gamma genus, and 22 representative HPV types for the species of Alpha and Beta genera. Full-lengths of E1-E2 ORF concatenations were also used in this study to explore possible incongruent tree topologies.
Models of nucleotide substitution were selected by the jModel Test version 3.7 [23], and the GTR+I+G substitution model was the best-fit model for both maximum likelihood and Bayesian analysis methods. ML phylogenetic trees were constructed by MEGA5 and data were bootstrap-resampled 1000 times to determine node support for the best trees. Bayesian phylogenetic analyses were performed with the MrBayes version 3.1.2 [24].The phylogenetic tree was displayed by using the Treeview [25].
PASC (Pairwise Sequence Comparison) is a web tool for analysis of pairwise identity distribution within viral families. The identities are pre-computed for every pair within the families with distribution plotted in the form of histogram where each bar corresponds to an interval of identities [26]. The global alignments function of this program was used to calculate the pairwise identities of these novel HPV types within the Papillomaviridae.

Ethics Statement
None of the 2210 samples from the Anyang population were collected solely for the purpose of this study. Samples of exfoliated skin cells were collected prospectively in our previous studies. All samples included in the study were collected in compliance with the Helsinki Declaration. This study was approved by the institutional review boards of Beijing Cancer Hospital. All participants gave written informed consent for use of their samples in research. In order to protect the rights of the participants, all samples used in the study were coded and tested anonymously.

Results
The clones and corresponding complete sequences were submitted to the International Reference Centre for

Genome organization of the eleven novel HPV types
The genome lengths of these 11 novel HPVs are 7128bp to 7415bp, and the genomic G+C content is between 35.83% (HPV 167) and 38.66% (HPV 166). They all have a genomic organization in common as has been reported in cutaneous PVs [10,12]. Each contains seven ORFs, which potentially encode five early genes E6, E7, E1, E2, E4, and two late genes L2 and L1 (Table 1). There are no nucleotide sequences for the E5 gene between early and late genes and such sequences are typically missing in β-PV and γ-PV [13]. The putative E4 ORFs of HPV 161, 162, 163, 165, 166, 167, and 170 do not have start codons, however sequences that may potentially encode this gene are present in these genomes, similar to E1 ORFs described previously [21]. Further validation of E4 transcripts should be carried out.
In silico analysis shows that the E6 ORF of all these novel HPV types contains two zinc-binding domains CxxC(x) 28-30 -CxxC, separated by 36 amino acids. The E7 ORF contains one zinc-binding domain [CxxC(x) 29-30 -CxxC], and six HPVs (HPV 161, HPV 162, HPV 165, HPV 166, and HPV 167) have the LxCxE motif for binding to the pRB protein, but the remainder of these novel HPV types lack this motif. The ATP binding site of the ATP-dependent helicase (GPPDTGKS) is conserved in the carboxy-terminal region of E1 in HPV 161, HPV 162, HPV 163, HPV 165, HPV 166, HPV 168, and KC5. The ATP binding site is modified slightly in HPV 164 and HPV 167 as "GPPDSGKS", in HPV 169 as "GVPDSGKS", and in HPV 170 as "GPPNTGKS" (the amino acid in bold indicate mutation site). The length of the upstream regulatory region (URR) is from 447 bp (HPV 165) to 581 bp (HPV 170) (Table 1), similar to the size of URRs in other HPV types of the genus Gammapapillomavirus [9,12]. Within the URR region of all these novel HPV types, putative E2-binding sites (ACCN 6 GGT) were identified (Table S1).

Phylogenetic representation of novel HPV types in relation to known Human Papillomavirus
Phylogenetic trees were produced using either the maximum-likelihood or Bayesian method. As very similar topologies were observed, we therefore combined the ML results and Bayesian results for the L1 tree and the E1-E2 tree, respectively ( Figure S1-S4). Phylogenetic analysis based on both the L1 gene and E1-E2 genes indicates that these 11 novel HPV types cluster into the Gamma genus, although there are topological inconsistences between early-gene-and lategene-derived phylogenies among several novel HPV types. Based on L1 ORF sequences, these 11 novel HPVs were classified into flowing species, as follows: HPV 163 to the γ-1 species; HPV 164 and HPV 168 to the γ-8 species; HPV 165 and KC5 to the γ-12 species; HPV 168 to the γ-4 species; HPV 169 to the γ-11 species; and HPV 170 to the γ-12 species. It is noteworthy that HPV 161, HPV 162, and HPV 166 may form a new monophyletic group, whose closest relationship is with the γ-15 species (Figure 2).
In the E1-E2 early-gene-derived phylogenetic analysis, HPV 161, HPV 162, HPV 166, and HPV 168 showed incongruent results compared with late-gene-derived L1 phylogenies. In the early gene analysis, HPV 161, HPV 162, and HPV 166 were monophyletic, but they clustered with the γ-12 species rather than the γ-15 species in the late gene analysis. In addition, the most closely related species of HPV 168 was γ-8 but not γ-4 based on the result of the late gene analysis.

Prevalence of novel HPV types in the healthy population
A total of 2210 skin exfoliated cell samples were collected from the Anyang rural population and tested for cutaneous HPV infection. The prevalence of novel HPV types was extremely low. The infection of novel HPV types appeared either alone or in combination with other HPV types (Table 3).

Discussion
The skin serves both as a barrier to prevent microbial infection and as a habitat for a diversity of commensal and pathogenic bacteria, fungus, and virus, which contribute to both human health and disease [27,28]. Human papillomavirus are distributed widely on human skin and exhibit great genetic diversity [6,8,29]. Unlike HPVs which inhabit the genital tract, HPVs from the cutaneous ecosystem engage in much milder relationships with human beings [4]. Present evidence suggests that infective cutaneotropic HPVs may only lead to epidermodysplasia verruciformis (EV) or other benign skin diseases [30][31][32]. Differences in the pathogenic capacity of HPVs are closely related to variations in the viral genome. The analysis of genomes provides precision sequence changes that may be evolutionarily important and comprise biologic forces driving HPVs into discrete ecosystems.
The majority of putative HPV types have been initially isolated from healthy skin without lesions suggesting the generally commensal nature of cutaneotropic HPVs [5,33,34]. The novel HPVs found in this study were all isolated from healthy individuals from a cohort study of esophageal cancer, in rural Anyang, Henan province, China [15]. Results of general HPV investigation in the Anyang rural population showed very low prevalence of these newly discovered HPVs. Generic primers rather than type-specific primers were used for detecting these novel HPVs so there is a possibility that the prevalence of these novel HPVs was underestimated.
Systematic positions of HPV 161 to HPV 170, and KC5 were defined by reconstruction of the phylogeny tree based on L1 in  68 HPV types. Although the L1 ORF is the gold standard for defining a PV type, there are many researchers who think that other genes also play important roles in classifying the phylogenetic relationships within the PV family [35,36]. Different genes within the same virus have different evolutionary rates, and different PVs have evolved at different rates [35,37]. In this study, the E1-E2 ORFs were employed in independent phylogenetic analysis to explore possible incongruent topologies. Both early-derived genes (E1-E2) and a late-derived gene (L1) demonstrated that all 11 novel HPVs belong to the γ-PV genus (Figure 1). However, there was incongruence in the topologies of four HPVs (HPV 161, 162, 166, and 168). Values for ML bootstrap support and Bayesian posterior probability were low in some branches inside the γ-PV, and this may be another possible reason for the incongruence between these two independent phylogeny analyses. The identification of phylogenetically representative HPV genome sequences is necessary to fill the gaps in our current knowledge about γ-PV genus.
In conclusion, eleven novel HPV types isolated from normal skin belong to the γ-PV genus. HPV 161, HPV 162, and HPV 166 are a monophyletic group and do not fit into any known species in γ-PV genus, suggesting that these may form a new HPV species in the γ-PV genus. The prevalence of these 11 novel HPVs in the normal population was very low, which strengthens the view that skin HPVs are highly diversified.