Genetic variation of human papillomavirus type 16 in individual clinical specimens revealed by deep sequencing

Abstract

Viral genetic diversity within infected cells or tissues, called viral quasispecies, has been mostly studied for RNA viruses, but has also been described among DNA viruses, including human papillomavirus type 16 (HPV16) present in cervical precancerous lesions. However, the extent of HPV genetic variation in cervical specimens, and its involvement in HPV-induced carcinogenesis, remains unclear. Here, we employ deep sequencing to comprehensively analyze genetic variation in the HPV16 genome isolated from individual clinical specimens. Through overlapping full-circle PCR, approximately 8-kb DNA fragments covering the whole HPV16 genome were amplified from HPV16-positive cervical exfoliated cells collected from patients with either low-grade squamous intraepithelial lesion (LSIL) or invasive cervical cancer (ICC). Deep sequencing of the amplified HPV16 DNA enabled de novo assembly of the full-length HPV16 genome sequence for each of 7 specimens (5 LSIL and 2 ICC samples). Subsequent alignment of read sequences to the assembled HPV16 sequence revealed that 2 LSILs and 1 ICC contained nucleotide variations within E6, E1 and the non-coding region between E5 and L2 with mutation frequencies of 0.60% to 5.42%. In transient replication assays, a novel E1 mutant found in ICC, E1 Q381E, showed reduced ability to support HPV16 origin-dependent replication. In addition, partially deleted E2 genes were detected in 1 LSIL sample in a mixed state with the intact E2 gene. Thus, the methods used in this study provide a fundamental framework for investigating the influence of HPV somatic genetic variation on cervical carcinogenesis.

Figure 1. Amplification of full-length HPV16 genomes by full-circle PCR.

(A) PCR was performed with PrimeSTAR^® GXL DNA polymerase and HPV16-specific primer-pairs as indicated. The amounts of HPV16/pUC19 used for the PCR template are also indicated above. The PCR products were analyzed by agarose gel electrophoresis. M, DNA size markers. (B) Scheme for full-circle PCR. PrimeSTAR^® GXL DNA polymerase generates short and long DNA products with primer-pair 1742F/1873R. (C) Full-circle PCR with DNA extracted from W12 cells, clone 20863 (high-copy HPV16 episomes) (lane 2) and clone 20850 (low-copy HPV16 episomes) (lane 3). M, DNA size marker (lanes 1) (D) Full-circle PCR using DNA isolated from 7 clinical specimens: 5 LSIL (lanes 2 to 6), and 2 ICC (lanes 7 and 8). M, DNA size marker (lanes 1).

Figure 2. Mutation frequency profile of full-length HPV16 genomes prepared from a plasmid.

The read sequences obtained with full-length HPV16 genomes prepared from a cloned plasmid were aligned to the reference HPV16 sequence, and mutation/error frequencies at each nucleotide position relative to the total coverage (number of reads that encompass each nucleotide, up to the maximum x 8,000) are presented in the landscape of the full-length HPV16 genome. A threshold line for a reliable mutation frequency (0.5%) is indicated with the red dotted line. The genome organization of HPV16 is indicated above: p97, the early promoter; p670, the late promoter; polyA(early) and polyA(late), the early and late polyadenylation signals, respectively.

Figure 3. Mutation frequency profile of full-length HPV16 genomes in W12 cells.

The read sequences obtained with full-length HPV16 genomes prepared from W12 cells were aligned to the reference HPV16 sequence (AF125673), and mutation/error frequencies at each nucleotide position are presented in the landscape of the full-length HPV16 genome. A threshold line for a reliable mutation frequency (0.5%) is indicated with the red dotted line. The genome organization of HPV16 is indicated above.

Figure 4. Mutation frequency profile of full-length HPV16 genomes in clinical specimens.

The read sequences obtained with full-length HPV16 genomes prepared from clinical specimens (2 LSIL samples, #2 and #3; 1 ICC sample, #6) were aligned to their de novo assembled complete genome sequences, and mutation/error frequencies at each nucleotide position are presented in the landscape of the full-length HPV16 genome. A threshold line for a reliable mutation frequency (0.5%) is indicated with the red dotted line. Peaks above 0.5% are indicated with red asterisks. The genome organization of HPV16 is indicated above.

Figure 5. Characterization of E1 variants by transient HPV16 replication assay.

(A) Left panel: schematic representation of the HPV16 E1 protein. The prototype E1 protein has leucine at position 294 and isoleucine at position 326, both of which are located in the DNA-binding domain, and glutamine at position 381 located in the oligomerization domain. Right panel: summary of aa variations found in the E1 protein from the clinical samples (#1 to #7). Proto, the E1 protein of European prototype. (B) E1-dependent replication of HPV16 origin-positive plasmid. The HPV16 origin-positive plasmid expressing Firefly luciferase (Fluc) and the origin-deficient plasmid expressing Renilla luciferase (Rluc) were transfected into HEK293 cells with increasing amounts of the prototype E1 expression plasmid together with the E2 expression plasmid. At 72 h after transfection, the ratio of the two luciferase activities (Fluc/Rluc) was measured as levels of the origin-dependent replication. The origin-deficient plasmid expressing Firefly luciferase (pGL4.50) was used as a negative control for HPV replication. Statistically significant differences (Welch’s t-test, p<0.01) are indicated with *. (C) Replication activity of E1 variants. Increasing amounts of expression plasmids for FLAG-tagged prototype or variant E1 proteins were transfected into HEK293 cells, and the levels of replication were measured at 72 h after transfection. Error bars represent the standard deviation of triplicate transfections. Statistically significant differences (Welch’s t-test, p<0.01) are indicated with *. The data are representative of three independent experiments. (D) Western blot analysis of E1 variants. FLAG-tagged prototype or variant E1 proteins expressed in HEK293 cells were detected with anti-FLAG antibody (Sigma-Aldrich) and the ECL prime Western blotting detection reagent (GE Healthcare, Buckinghamshire, England). Tubulin was visualized with anti-α/β-tubulin antibody (Cell Signaling, Danvers, MA) as loading control.

Figure 6. HPV16 genome deletions in a clinical specimen revealed by read mapping.

Mapping of paired-end read sequences to de novo assembled full-length HPV16 genome sequence for LSIL sample 2 is visualized by Tablet. Deletions in the E2/E4 region (A) and in the region between E5 and L2 (B) are shown. (C) The genome organization of HPV16. The deleted regions in E2/E4 and between E5 and L2 are indicated with red and blue dotted-line boxes, respectively.