Viral Metagenomics Reveals a Putative Novel HPV Type in Anogenital Wart Tissues

Viral metagenomics is widely employed to identify novel viruses in biological samples. Recently, although numerous novel human papillomavirus (HPV) types have been identified in clinical samples including anogenital warts (AGWs), many novel HPV sequences remain to be discovered. In this study, a putative novel HPV type designated as HPV-JDFY01 was discovered from library GW05 with 63 sequence reads by the viral metagenomic technique. Its complete genomic sequence was determined by PCR to bridge the gaps between contigs combining Sanger sequencing. The complete genome of HPV-JDFY01 is a 7186 bp encoding 7 open reading frames (ORFs) (E6, E7, E1, E2, E4, L2 and L1) and contains a 487 bp long control region (LCR) between L1 and E6. Sequence and phylogeny analysis indicated that HPV-JDFY01 shared the highest sequence identity of 74.2% with HPV-mSK_244 (MH777383) and well clustered into the genus Gammapapillomavirus. It has the classical genomic organization of Gammapapillomaviruses. Epidemiological investigation showed that one out of the 413 AGW tissue samples was positive for HPV-JDFY01. Further research with large size and different type of samples should be performed to elucidate the epidemiologic status of HPV-JDFY01.


Introduction
HPVs are present as small, nonenveloped and circular double-stranded DNA viruses with an icosahedral capsid, and they consist of an average genome length of 8000 bp [1,2]. The whole genome can be broadly grouped into early regions, late regions and a long control region according to their different functions [3][4][5]. The outer capsid of HPV is comprised predominantly of seventy-two major capsid protein L1, which plays an essential role in the attachment during the viral entry and the assembling of virions [6,7]. L1 protein has also been applied for the classification or defining of a novel HPV type if the variability proportion is over 10% [3].
Gammapapillomavirus accounts for over 50% of total HPVs and owns 83.3% of the novel types from 2014, which demands significant attention [8]. Little is known about its genetic characteristics, but it has been reported that all of the Gamma-6 lack the E6 protein to promote carcinogenesis [9]. In addition, patients with warts, hypogammaglobulinemia, infections or myelokathexis (WHIM) syndrome seem to be inclined to be infected by Gammapapillomavirus [10]. In the past four years, only four novel HPV types had been identified, including HPV226 (Gamma-6), HPV227 (Beta-2), HPV228 (Gamma-27) and HPV229 (Gamma-7) (https://www.hpvcenter.se/), three of which belonged to Gammapapillomavirus. Latsuzbaia et al. [11] discovered that the strain HPV226 contained one TATA box and three palindromic E2-binding sites in LCR, one zinc-binding domain in E7 and one ATP-binding site in E1. Other special structural features, such as polyadenylation sites and LxCxE motif were identified in HPV227 [12]. Both HPV226 and HPV227 were detected by next-generation sequencing (NGS) after rolling circular amplification [11,12]. Viral Metagenomics is based on NGS and can detect all the viral nucleic acid directly in biological samples so as to effectively find possible new viruses. In our study, using viral metagenomic method, a putative novel HPV (HPV-JDFY01) was discovered from AGWs tissue samples and the complete genome was acquired by combining PCR and Sanger sequencing. Sequence and phylogeny analysis were further performed to fully characterized this putative novel type of HPV.

Sample Collection and Preparation
A total of 110 tissue samples of wart from AGW patients, including 34 females and 76 males, were collected from September 2016 to August 2017 from the Department of Dermatology and Venerology at the Affiliated Hospital of Jiangsu University. The mean age of AGW patients was 42.5 years and the mean disease duration was 3.26 months. Each sample was transferred into 1.5 mL test tube with 1 mL dulbecco's phosphate-buffered saline (DPBS) and then added into several 5 mm sterile stainless-steel beads to fully grind tissue for 5 min. After freeze and thaw 3 times on dry ice for 3 min each time, the tissue supernatants were obtained after centrifuging for 10 min at 12,000× g under the conditions of 4 • C.

Viral Metagenomic and Bioinformatics Analysis
The 110 supernatants of wart tissues were randomly pooled into 11 sample pools, each containing 10 samples. These sample pools were filtered through a 0.45 µm filter (Merck Millipore Ltd., Tullagreen, Carrigtwohill, Darmstadt, Germany, Co. Cork, IRL) to remove cellular particles in the sample [13,14]. The filtrates were digested by DNase and RNase at 37 • C for 60 min to remove non-viral encapsidated nucleic acids [15,16]. Then, a QIAamp Viral RNA Mini Kit (QIAGEN GmbH, QIAGEN Strasse 1, 40724 Hilden, Germany) was utilized to extract the remaining nucleic acids in the samples [17]. Reverse transcription (RT) and Klenow reactions were conducted to synthesize double-stranded DNA, which was prepared for constructing the 11 cDNA libraries by Nextera XT DNA Sample Preparation Kit (Illumina, Inc., San Diego, CA 92122 USA) [18,19]. Libraries were subjected to NGS on the MiSeq Illumina sequencing platform. The generated NGS data was processed by in-house analysis pipeline. Briefly, the bioinformatics steps were performed on a cluster of 32 nodes running Linux, and taking the Phred Score Q30 as the threshold trimmed the low sequencing quality tails [20]. Adaptors and primer sequences were trimmed by VecScreen in the default parameters [21]. These cleaned reads from sequencing were de novo assembled into contigs within each barcode group by Ensemble assembler [22]. The assembled contigs, along with unassembled reads, were matched against viral proteome database by BLASTx with an E-value cutoff of <10 −5 , where the viral proteome databased included NCBI virus reference proteome together with viral proteins sequences from NCBI nr fasta file.

Acquiring of Complete Genome and PCR Reactions
Finding six contigs ranging from 123 bp to 645 bp in length in GW05 pool, which showed low sequence similarity to known HPVs, 2 sets of specific PCR primers derived from a contig belonging to L1 gene of HPV were, respectively designed to confirm the existence of the putative novel HPV sequence in samples from GW05 pool. PCR primer  Table 1. The screening conditions of PCR were as follows: 95 • C for 5 min, 31 cycles of 95 • C for 30 s, 48 • C (for the first round) or 52 • C (for the second round) for 30 s and 72 • C for 40 s and a final extension at 72 • C for 5 min. The circular complete genome of HPV-JDFY01 was acquired by PCR primers bridging the gaps between contigs and the PCR products were cloned and subjected to Sanger sequencing. In order to assess the prevalence of HPV-JDFY01 in AGW samples, an epidemiological investigation was performed by PCR screening method described above for a total of 413 AGW specimens collected from August 2017 to August 2020.

Phylogenetic Analysis
Phylogenetic analysis was performed based on amino acid sequences of predicted ORFs to analyze the genetic relationships between HPV-JDFY01 and the other related HPV types which included the representative types of Alphapapillomavirus, Betapapillomavirus and Gammapapillomavirus. Multiple alignments were conducted using ClustalW in MEGA7 with the default settings [23]. The aligned file was converted to Nexus (.nex) formmat by Geneious 11.1.2 and subsequently modified by Notapade++ [24]. The mean standard deviation of split frequencies was controlled as low as possible when generating phylogenetic trees by Mrbayes 3.2.7 with Bayesian inference (BI) [25].

Results
The NGS results indicated that the 11 libraries generated a total of 6,509,970 reads, 726,777 (11.16%) of which showed significant sequence similarity to HPV. The library GW02 generated the largest number of total reads (1,300,146) while the library GW09 contained the most abundant HPV-related sequence reads (227,374). Among the HPV sequence reads in these libraries, Alpha-10 (26,708) ranking first was followed by Alpha-8 (20,929), Alpha-14 (1622) and Alpha-9 (1066). HPV11 occupied the dominant part in Alpha-10 (95.41%) and Pathogens 2022, 11, 1452 4 of 10 HPV of the GW05 library (33.87%), respectively, in accordance with that of the literature reported previously [26,27]. HPV6, another common subtype of AGW, accounted for 0.42% in GW05 library. Except for unclassified papillomaviridae, the other HPV types identified in the 110 wart samples were presented in Table S1. Given that six contigs ranging from 123 bp to 645 bp in length in the GW05 pool showed low sequence similarity to known HPVs, 10 individual samples in this pool were screened by the nested PCR. The result indicated only one sample was positive. PCR combining Sanger sequencing was then used to bridge the sequence gaps and obtain the whole genome of the putative novel HPV eventually.

Complete Genomic Structure of HPV-JDFY01
A putative novel HPV type (HPV-JDFY01) was discovered from AGW, of which the complete circular genome is 7186 bp in length with a GC content of 36.6%. The details of the seven predicted ORFs are described as follows: Early proteins: E6 (nucleotide

Phylogenetic Analysis of HPV-JDFY01
The amino acid sequence of HPV-JDFY01 L1 was aligned with 11 typical reference sequences from Alphapapillomavirus, 9 from Betapapillomavirus, 15 from Gammapapillomavirus and 13 from unclassified papillomavirus so as to generate the phylogenetic tree that clearly reflected the evolutionary relationships of these representative HPV types. From the results of the phylogenetic analysis, the HPV-JDFY01, together with HPV-mSK_244 ( Figure 2, Figures S1-S5) was clustered into the same clade of Gammapapillomavirus, a genus with the highest number of novel HPVs in recent years [8]. Both HPV-JDFY01 and HPV-mSK_244 were directly adjacent to the HPV112 which belongs to the Gamma-8 species. The sequence identity and phylogenetic analysis supported that HPV-JDFY01 was closely related to Gamma-8 types.

Phylogenetic Analysis of HPV-JDFY01
The amino acid sequence of HPV-JDFY01 L1 was aligned with 11 typical reference sequences from Alphapapillomavirus, 9 from Betapapillomavirus, 15 from Gammapapillomavirus and 13 from unclassified papillomavirus so as to generate the phylogenetic tree that clearly reflected the evolutionary relationships of these representative HPV types. From the results of the phylogenetic analysis, the HPV-JDFY01, together with HPV-mSK_244 (Figure 2, Figures S1-S5) was clustered into the same clade of Gammapapillomavirus, a genus with the highest number of novel HPVs in recent years [8]. Both HPV-JDFY01 and HPV-mSK_244 were directly adjacent to the HPV112 which belongs to the Gamma-8 species. The sequence identity and phylogenetic analysis supported that HPV-JDFY01 was closely related to Gamma-8 types.

Epidemiological Analysis of HPV-JDFY01
Epidemiological investigation showed that only one out of the 413 AGW tissue samples was positive for HPV-JDFY01. The typically clinical and pathological manifestations were shown in Figure 3A,B.

Epidemiological Analysis of HPV-JDFY01
Epidemiological investigation showed that only one out of the 413 AGW tissue samples was positive for HPV-JDFY01. The typically clinical and pathological manifestations were shown in Figure 3A,B.

GenBank Accession Number
The complete genomic sequence of HPV-JDFY01 was submitted to the GenBank database under accession number MW311489.

Discussion
Viral metagenomics, based on NGS, can detect nucleic acid of all viral pathogens from samples and have the potential to become crucial for clinical diagnosis in the near future [35][36][37][38]. We report a genome of HPV-JDFY01, a putative novel type in Gammapapillomavirus, identified by viral metagenomics from AGWs collected in Jiangsu Province, China.
Previous literatures considered that HPVs of the same species share ≥ 71% L1 similarity [1,3]. Phylogenetic analysis indicated that HPV-JDFY01 shared the highest sequence identity of 74.2% with HPV-mSK_244 based on the nucleotide sequence of L1. Thus, HPV-JDFY01 could be assumed as a putative novel HPV type and may have similar biological functions with HPV-mSK_244. HPV-mSK_244 was identified from 207 skin and nares swabs from DOCK8-deficient patients, who have recurrent cutaneous and systemic infections and can suffer from the predisposition and susceptibility to cancer [7]. However, HPV-JDFY01 was found in immunocompetent AGW patients. Additionally, in the phylogenetic analyses of L1 protein, HPV-JDFY01 and HPV-mSK_244 were clustered together with HPV112, which belongs to the Gamma-8 species. It is worth noting that HPV112, a member of the Gamma-8 species, was also cloned from a 25-year-old male with condyloma acuminata [39]. The other two types of Gamma-8, HPV164 and HPV168, were isolated from exfoliated skin cell samples of healthy individuals [40,41]. In summary, the above evidences indicate that HPV types from the Gamma-8 species may engage in much milder relationships with human beings. Additionally, the relative abundances of DNA viruses, mostly Papillomaviridae and Polyomaviridae, were higher on the skin than in the oral cavity and feces of DOCK8-deficient patients [7]. HPV-JDFY01 was isolated from tissue samples of wart while HPV112, HPV164 and HPV168 were all isolated from skin [39,40]. Therefore, it can be stated that the members of Gamma-8 tend to infect skin and may be distinguished into cutaneotropic HPV types.
Even though further studies are needed to assess the relationship between HPV-JDFY01 and Gamma-8, the genetic characteristics of HPV-JDFY01 were further analyzed by comparing with the previous literatures, which were related to other HPV members of Gamma-8. We found that most of specific characteristics of Gamma-8 could be identified in HPV-JDFY01, including zinc-finger domains, PDZ binding domains, ATP-binding sites, nuclear localization signals, nuclear export signals, NLS-like signals, furin cleavage

GenBank Accession Number
The complete genomic sequence of HPV-JDFY01 was submitted to the GenBank database under accession number MW311489.

Discussion
Viral metagenomics, based on NGS, can detect nucleic acid of all viral pathogens from samples and have the potential to become crucial for clinical diagnosis in the near future [35][36][37][38]. We report a genome of HPV-JDFY01, a putative novel type in Gammapapillomavirus, identified by viral metagenomics from AGWs collected in Jiangsu Province, China.
Previous literatures considered that HPVs of the same species share ≥ 71% L1 similarity [1,3]. Phylogenetic analysis indicated that HPV-JDFY01 shared the highest sequence identity of 74.2% with HPV-mSK_244 based on the nucleotide sequence of L1. Thus, HPV-JDFY01 could be assumed as a putative novel HPV type and may have similar biological functions with HPV-mSK_244. HPV-mSK_244 was identified from 207 skin and nares swabs from DOCK8-deficient patients, who have recurrent cutaneous and systemic infections and can suffer from the predisposition and susceptibility to cancer [7]. However, HPV-JDFY01 was found in immunocompetent AGW patients. Additionally, in the phylogenetic analyses of L1 protein, HPV-JDFY01 and HPV-mSK_244 were clustered together with HPV112, which belongs to the Gamma-8 species. It is worth noting that HPV112, a member of the Gamma-8 species, was also cloned from a 25-year-old male with condyloma acuminata [39]. The other two types of Gamma-8, HPV164 and HPV168, were isolated from exfoliated skin cell samples of healthy individuals [40,41]. In summary, the above evidences indicate that HPV types from the Gamma-8 species may engage in much milder relationships with human beings. Additionally, the relative abundances of DNA viruses, mostly Papillomaviridae and Polyomaviridae, were higher on the skin than in the oral cavity and feces of DOCK8-deficient patients [7]. HPV-JDFY01 was isolated from tissue samples of wart while HPV112, HPV164 and HPV168 were all isolated from skin [39,40]. Therefore, it can be stated that the members of Gamma-8 tend to infect skin and may be distinguished into cutaneotropic HPV types.
Even though further studies are needed to assess the relationship between HPV-JDFY01 and Gamma-8, the genetic characteristics of HPV-JDFY01 were further analyzed by comparing with the previous literatures, which were related to other HPV members of Gamma-8. We found that most of specific characteristics of Gamma-8 could be identified in HPV-JDFY01, including zinc-finger domains, PDZ binding domains, ATP-binding sites, nuclear localization signals, nuclear export signals, NLS-like signals, furin cleavage motifs, polyadenilation sites, palindromic E2-binding sites and TATA boxes. In previous literature, it has been reported that the E7 protein of HPV112 contained a pRB protein binding domain, although it was not totally conserved [39]. However, pRB protein binding motif (LxCxE) have not been detected in E7 protein of HPV-JDFY01. This means lacking the LxCxE motif may have the effect of attenuating the association between E7 and pRB [42] and indicate transforming properties of HPV-JDFY01. A highly conserved NLS motif can promote nuclear localization and bind to the nuclear matrix [43], which was also not identified in E2 of HPV-JDFY01. Moreover, three internal PDZ-binding motifs (PTAL, LSLV and ESLV) were identified in the putative E6 protein. Nevertheless, these motifs may not possess biological functions [44]. Further structural analysis will be needed to better understand the relationship between HPV-JDFY01 and the other members from Gamma-8.
On the basis of the highly conserved, immunogenic and self-assembly natures in HPV L1 protein, virus-like particles (VLPs) of HPV have been designed and developed for the purpose of prophylactic HPV vaccines [45]. Further prediction analyses of advanced structure for novel HPV L1 protein are thus warranted.

Conclusions
In this study, a putative novel Gammapapillomavirus named HPV-JDFY01 was detected by a viral metagenomics method from AGWs, whose complete genome was determined and genomic structure and phylogeny were fully characterized. The functional, biophysical and carcinogenic properties of HPV-JDFY01 need to be elucidated in the future research.

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki, and approved by Medical Ethical Committee at the Affiliated Hospital of Jiangsu University (Reference code: ujsah2016219).

Informed Consent Statement:
Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.