Phylogeny and Classification of Human Papillomavirus (HPV)16 and HPV18 Variants Based on E6 and L1 genes in Tunisian Women with Cervical Lesions

Background: High-risk human papillomavirus (HPV) types are the main etiological factors for cervical cancer.
HPV16 and HPV18 are generally the most common forms associated with development of high-grade cervical lesions.
This study was undertaken to identify intratypic variants of HPV16 and HPV18 among women with cervical lesions
in Tunisia. Materials and Methods: DNA was extracted from cervical samples collected from 49 women. using a
PureLinkTM Genomic DNA mini Kit (Invitrogen). E6 and L1 open reading frames (ORF) were amplified by PCR
and viral DNA amplicons were subjected to automated sequencing using Big Dye Terminators technology (Applied
Biosystems). The obtained sequences were analyzed using an appropriate software program to allow phylogenetic
trees to be generated. Results: HPV16 and HPV18 were detected in 15 and 5 cases, respectively. HPV16 E6 sequences
clustered with the European German lineage (A2) whereas one isolate diverged differently in the L1 region and
clustered with the African sub-lineage (B1). HPV 18 E6 sequences clustered with the European sub-lineage (A1)
but L1 sequences clustered as a new clade which diverged from A1-A5. Conclusions: Our results suggest that the
distribution of HPV16 and HPV18 sequences in women with cervical lesions in Tunisia is mainly related to European
epidemiological conditions and point to the presence of recombinant HPV forms.


Phylogeny and Classification of Human Papillomavirus (HPV)16 and HPV18 Variants Based on E6 and L1 genes in
The prevalence of genital HPV infection in women differs considerably between countries and regions (Clifford et al., 2005;Seoud, 2012) as well as between different risk groups. In Tunisia, a 13.2% prevalence in the Grand Tunis region has been reported (Ardhaoui et al., 2016). In studies conducted on prostitutes, HPV-DNA prevalence was 39%-45% (Hassen et al., 2003;Demarco et al., 2006). HPV typing in these high-risk groups as well as in Tunisian women with precancerous and cancerous lesions showed the predominance of HPV16 as expected from the data of other countries (Missaoui et al., 2010).
It was reported that HPV16 and HPV18 variants are reflective of the migration patterns of Homo sapiens and suggested that HPV16 and HPV18 variant lineages might have diverged through genetic isolation at approximately the same time as Homo sapiens began establishing residence in different continental regions (Fu Xi et al., 2006;Chen et al., 2005;Chen et al., 2009;Pimenoff et al., 2016). Indeed, HPV16 and HPV18 isolates worldwide can be classified into 4 and 3 lineages, respectively, according to nucleotide variations in their genomes (Chen et al., 2009). To our knowledge, only two reports on intratypic HPV variants in Tunisia have been published to date and HPV16 and HPV18 variants of European sublineage were found to be the most prevalent (KrennHrubec et al., 2011;Ghedira et al., 2016). It has been reported that HPV intratypic variants differ in oncogenic potential (Villa et al., 200;Hang et al., 2016). In the present study in order to identify the HPV16 and HPV18 intratypic variants in women with cervical disease, we have sequenced and analyzed the E6 and L1 open reading frames (ORF) of both viruses identified among a sample of Tunisian women with precancerous or cancerous cervical lesions.

Patients and cervical samples
The study population consisted of 49 Tunisian women with cervical lesions consulting either at the Department of Gynecology and Obstetrics, La Rabta Maternity and Neonatology Center, Tunis, or in a private gynecologic office. Their mean age was 45 year (range, 21-64 years). A questionnaire about personal information was filled for each woman; clinical characteristics of patients are summarized in Table 1. Cervical samples were collected in SurePath™ vials (BD Diagnostics, France) for Pap smear assay and HPV typing. The study was conducted in compliance with the Helsinki Declaration and approved by the Habib Thameur Hospital local Ethical Committee, Tunis. Tunisia.

DNA Extraction and HPV amplification
DNA was extracted from cervical samples using PureLinkTM Genomic DNA mini Kit (Invitrogen, Carlsbad, CA, USA) according to manufacturer's instructions. DNA extracts were stored at -20°C until use.
HPV-DNA amplification was performed in the presence of consensus primers for genital HPVs and in the presence of type-specific primers for HPV 16 and HPV 18. The primer pair MY09/MY11 (Manos et al., 1989) was used as external primers followed by GP5+/GP6+ (Snijders et al., 1990) as internal primers for nested PCR to amplify a region of about 150 bp in ORF L1. Type-specific PCR was performed in the presence of either primer pair E6-HPV16U/E6-HPV16D for HPV 16 or E6-HPV18U/ E6-HPV18D for HPV 18. These primers specify a region of 420 bp for HPV16 or 456 bp for HPV18. Primer sequences are shown in Table 2.
Target DNA was amplified in a 50 µl reaction mixture containing: 1x Taq polymerase buffer, 1.5 mM MgCl2, 200µM each dNTP, 1.25 units HotStart Taq DNA polymerase (Invitrogen, Courtaboeuf, France), 0.5 µM each primer, 10µl DNA extract for simple PCR or 2 µl of first round PCR product for nested PCR. Thermocycler programs were as follows. For simple PCR with MY09/ MY11 and for type-specific PCR: 40 cycles of denaturation at 94°C for 1 minute, annealing for 1 minute at 55°C, extension at 72°C for 1 minute and a final extension step at 72°C for 7 minutes. For the second round of nested PCR the temperature program was as above except for annealing which was performed at 50°C. PCR reactions were carried out in a GeneAmp PCR system 9700 Thermal Cycler (Applied Biosystems). Each PCR run included a blank control containing water instead of target DNA, a negative control DNA from HEp-2 cells and a positive control containing either HPV 16 DNA from SiHa cells or HPV18 DNA from HeLa cells (Adey et al., 2013). The quality of DNA extracts was assessed by amplification of a 536 bp sequence of the β globin gene using the primer pair KM29/RS42. All PCR products were visualized on 2% agarose gels in the presence of ethidium bromide.

HPV Sequencing and Phylogenetic Analysis
Sequencing of the purified PCR products was carried out in the forward and reverse directions using the

Results
Among the 49 woman included in the study, 15 (30.6%) and 5 (10.2%) were positive for HPV16 and HPV18, respectively and the remaining 29 samples were positive for other HPV types. As shown in Table 4, 65% (13/20) of the women positive for HPV16 or HPV18 had precancerous (HSIL) or cancerous lesions.
All PCR products were sequenced and sequences from each region were submitted to GenBank under accession numbers KT931993 to KT932004, KU175625 to KU175627 and KX75949to KX759668.
Variant distribution was determined through L1 and E6 ORF sequences. HPV16 isolates were initially examined for nucleotide variation within a 114-nucleotide (nt) sequence in ORF L1 (nt 6654 to 6767). Sequence polymorphism was identified at 1 nt position through the ORF L1. All HPV16 strains appeared A2-like in E6 but one strain (MJR9) diverged differently within the L1 region. Indeed, within the 392-nt sequence in ORF E6 (nt 228 to 620), HPV16 strains showed a T350G nt change which is specific to the European German lineage A2 [27] and lead to a L83V substitution, whereas in ORF L1, one corresponding primer pairs and Big Dye Terminators Cycle Sequencing Ready Reaction kit (Applied Biosystems, Foster City, CA, USA). Products of the sequencing PCRs were analyzed in an automatic ABI Prism 3100 Avant-Genetic Analyser (Applied Biosystems).
The obtained sequences were analyzed, using BioEdit (Hall, 1999) and Mega 6 (Tamura et al., 2013) softwares. Multiple alignments were generated by ClustalW integrated into BioEdit. Phylogenetic trees were generated using TreeView and DNA Maximum Likelihood program with molecular clock (DNAMLK) in Mega 6. The obtained isolates were aligned against references sequences shown in Table 3.   isolate showed the (G6719A) substitution, only found in the African sub-lineage B1. All HPV18 strains were similar to the European sublineage (A1-A5) through the E6 region, but interestingly, all the HPV18 isolates showed a G to A silent transition through ORF L1 at position 6731.

Phylogenetic analysis
Four phylogenetic trees based on L1 and E6 nt sequences for HPV16 and HPV18 were constructed using Maximum Likelihood Tree integrated into MEGA6 program (Figure 1).
For HPV16, the phylogenetic tree corresponding to E6 (Figure 1a) showed that all isolates clustered with European German sub-lineage (A2), whereas the phylogenetic tree corresponding to L1 HPV16 (Figure 1b) showed that fourteen HPV16 isolates clustered closely with European variants (A1-A3), while one isolate (MJR9) clustered with African sub-lineage (B1). This suggests that MJR9 could be a recombinant between European and African variants.
For HPV18, the phylogenetic tree corresponding to E6 (Figure 1c) showed that all isolates clustered with the European variant (A1), while for L1, phylogenetic tree ( Figure 1d) showed that these five isolates constitute a new clade due to G to A silent transition at position 6731, suggesting the possibility of a hitherto undescribed HPV18 variant.

Discussion
Intratypic variants of HPV16 and HPV18 show a particular geographical and ethnical distribution worldwide (Ong et al .,1993;Cornet et al., 2012 ;Chen Figure 1. Phylogenetic Trees of HPV16 and HPV18 Based on L1 and E6 Nucleotide were Constructed Using Maximum Likelihood Tree Integrated into MEGA6 Program. et al., 2015).
Contrary to the previous studies from Tunisia (Hassen et al., 2003;Ghedira et al., 2016) women enrolled in the present study were selected on the basis of cervical lesions. Therefore, the high frequency of precancerous or cancerous cervical lesions in these women positive for HPV16 or HPV18 is not surprising and, as expected, HPV16 was the most frequent type. Ennaifer et al., (2011) have reported predominance of HPV16 and HPV18 (92.1%) in cervical cancer cases from Tunisian women (Ennaifer et al., 2011). A similar result has been reported from a neighbor country of Algeria with a prevalence of 61.4% for HPV16 and 15.8% for HPV18 (Hammouda et al., 2011).
All HPV16 and HPV18 isolates identified in this study were closely related to European lineages. This finding tallies with worldwide reported data (Li et al., 2011;Bravo and Felez-Sanchez., 2015;Kabekkodu et al., 2015;Pimenoff et al., 2016, KrennHrubec et al., 2011. The dominance of HPV16 European variants among this group of Tunisian women is in agreement with data obtained in Morocco showing the predominance (58.3%) of European variant among Moroccan woman (Qmichou et al., 2013). However the proportion of African variants (31.1%) was higher in Morocco than in the present and a previous Tunisian study (Ghedira et al., 2016). However, the observation of a probable European/African HPV16 recombinant might reflect the contribution of populations from sub-Saharan African origin in the epidemiology of HPV in Tunisia. On the other hand, the fact that all HPV18 isolates constituted a new clade suggests a particular virus host co-evolution in this North African country.
It has been reported that intratypic HPV variants present a different oncogenic potential (Villa et al., 2000;Hang et al., 2016). In the present study, the homogeneity in HPV strains precluded to perform this analysis. Further studies will be necessary to determine whether the HPV16 recombinant or the HPV18 variants identified in this study have a particular oncogenic potential.
This study is the first to date to propose an evolutionary scenario for both HPV16 and HPV18 in Tunisia. The current study provides a phylogenetic classification among a group of Tunisian woman and confirms that HPV16 and HPV18 infections in Tunisia are much more due to European variants than to African variants. This suggests that the distribution of HPV16 and HPV18 are closely related to the European pattern and could reflect historical links and migrations between Tunisia and the European continent.
Further studies are required to analyze the relationship between HPV16 and HPV18 Tunisian variants and viral pathogenicity and the possible impact of the variability observed in ORF L1 on the HPV vaccine efficacy.