Rearrangements of viral and human genomes at human papillomavirus integration events and their allele-specific impacts on cancer genome regulation

Long-read whole-genome sequencing of cervical tumors

We sequenced 72 cervical cancer samples from the HTMCP cohort using whole-genome Oxford Nanopore Technologies (ONT) long-read sequencing (Methods). The samples yielded an average of 102 gigabase pairs (Gb) of data (range 52.6–153 Gb; median coverage = 34×; Supplemental Fig. S1). We achieved long read lengths in our samples (median N50 = 17.5 kb; range 9.0–34.1 kb; Supplemental Fig. S1) that were of high read quality, as measured by the base error rate (median = 4.8%; range 2.2%–8.4%) and the proportion of artifactual chimeric DNA fusions (median = 5.0%; range 1.0%–12.4%). These data are supplemented by high-quality, short-read whole-genome and RNA sequencing (RNA-seq), as previously described (Gagliardi et al. 2020).

Detection of HPV integration using long-read sequencing

The molecular and clinical properties of the 72 samples are shown in Figure 1A (Supplemental Table S1). From a sequence perspective, we define HPV integration as the recombination of HPV with human DNA, which can be detected from chimeric reads that align to both the HPV and human genomes. We observed that HPV integration events involved at least two double-stranded human genome breaks and two double-stranded HPV genome breaks (Fig. 1B). Breakpoints were grouped into an “integration event” if either of the following conditions was met: (1) the HPV breakpoints co-occurred on one or more of the same reads, or (2) the breakpoints mapped within 500 kb of each other (The Cancer Genome Atlas Research Network 2017; Gagliardi et al. 2020; Methods). Applying these conditions, we detected 438 integration breakpoints belonging to 129 integration events across 69 samples with HPV integration (Fig. 1A; Supplemental Table S2). To determine whether an event was transcribed, we searched for HPV–human fusion transcripts occurring in the vicinity of the integration event using existing RNA-seq data (Methods). At most, two events per sample produced fusion transcripts, but even samples without a recurrent RNA fusion site near an integration event expressed HPV (Fig. 1A; Supplemental Table S2). In three samples with deep and evenly distributed read coverage across the HPV genome, we did not detect HPV integration breakpoints, indicative of highly amplified episomal HPV. One sample harbored three integration events from two different HPV types, HPV16 (two events) and HPV59 (one event), although only the HPV59 event was expressed (Supplemental Table S2).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 1.

Detection and categorization of HPV integration events in cervical cancers. (A) The number of HPV integration events detected in the DNA (# of events) and detected in DNA and RNA (# of expressed events) across the HTMCP samples, as well as clinical and molecular characteristics. (B) Schematic illustrations of integration event categories. Orange segments indicate the insertion of HPV integrants. (C) Decision chart for categorizing HPV integration events (“bps” = breakpoints). (D) The frequency of the HPV integration categories across the samples. The percentage of events that produce an HPV–human fusion transcript is indicated for each integration type. (E) The percentage of events belonging to each integration category for HPV16, HPV18, HPV45, and all other HPV types. (F) The genomic locations of integration events across the cohort, colored by the transcriptional status of the event. Bins with 2+ integration events are highlighted with boxes. Notable cancer genes within bins recurrently affected by integration are indicated. (G) Gene expression differences between samples with HPV integration within 1 Mb of MYC, TP63, and KLF5 compared to the remaining samples in the data set. Box plots represent the median and upper and lower quartiles of the distribution; whiskers represent the limits of the distribution (1.5 IQR below Q1 or 1.5 IQR above Q3). P-values were calculated using the Wilcoxon rank-sum test.

Next, we compared the number and genomic locations of HPV integration breakpoints identified in our long-read sequencing data to previous short-read sequencing data (Gagliardi et al. 2020; Methods). Across the 69 HTMCP samples with HPV integration, 424 and 438 integration breakpoints were called in the short- and long-read sequencing data sets, respectively (Supplemental Tables S2 and S3). Integration breakpoint calls using short- and long reads were identical for 33 samples. Eighteen samples contained more calls in the long-read data, and 18 contained more calls in the short-read data (Supplemental Fig. S2A,B). Two samples where HPV integration was not detected with short reads contained one integration breakpoint each in the long-read data. One of these samples had even and deep read coverage indicative of episomal HPV, but an integration event was detected at a low variant allele frequency (VAF = 0.00025), suggesting that integration may have occurred subclonally (Supplemental Table S2). The second sample contained an integration site within a repetitive region of Chromosome 21p11 that was detected with long reads, but the region contained no short reads with adequate mapping quality.

Two samples had highly discordant integration breakpoint calls when detecting HPV–human junction sites, with an HPV58-integrated sample having 59 more calls in the long-read data (82 vs. 23) and an HPV16-integrated sample having 53 more calls in the short-read data (56 vs. 3; Supplemental Fig. S2A,B). We confirmed the overlap between integration breakpoint calls in each technology for these two samples. In the HPV58-integrated sample, 67/105 (64%) integration breakpoints called in either the long- or short-read data overlapped repetitive regions, and 30/59 (51%) breakpoints resided in repetitive regions in the HPV16-integrated sample, suggesting that read alignments may be confounded by repetitive sequences.

As expected, the number of integration breakpoint calls per sample was correlated between technologies (Supplemental Fig. S2C; Spearman's correlation, R = 0.78, P = 3.8 × 10⁻¹⁵). Both technologies detected similar proportions of breakpoint calls within CpG islands and exonic, intronic, and repetitive regions (Supplemental Fig. S2D). Long reads detected more integration breakpoints per sample, and discordance between calls made using long- and short-read technologies could be due to low-confidence mapping of short reads in repetitive regions.

Structural alterations associated with HPV integration events

Based on the sequence alignment patterns of HPV-containing reads and the HPV integration breakpoints, we defined five categories of integration structures in the 69 integrated samples. These were (del)etion-like, (dup)lication-like, translocation, multi-breakpoint, and repeat region integrations, which we categorized using a custom workflow (Fig. 1B,C; Methods). Multi-breakpoint integrations were the most prevalent, comprising 41/129 (32%) of integration events in our data set (Fig. 1D; Supplemental Table S2). Next were dup-like (37/129) and del-like integrations (29/129), which differed by the presence or absence of the intervening human sequence between the two HPV breakpoints (Fig. 1B,C). We also investigated the potential of integration events existing as extrachromosomal circular DNA (ecDNA) using the events’ de novo assemblies (Methods). We identified eight events as potential ecDNAs, most of which (6/8) were categorized as dup-like integrations (Supplemental Table S4). Translocation (3/129), and repeat region (3/129) integrations were rare (Fig. 1D; Supplemental Table S2). The remainder of the events (16/129) were unmatched single breakpoints that we were unable to categorize (Supplemental Table S2). While multi-breakpoint events were the most prevalent, only 49% of these events produced HPV fusion transcripts, in contrast to 76% of dup-like events. In fact, dup-like events were more often expressed when compared to all other types of events (Fisher's exact test, P = 0.00096, Fig. 1D). The distributions of the integration categories also differed between HPV types (Fig. 1E).

Consistent with previous reports, we observed regions of recurrent HPV integration (>2 events within 1 Mb bins) in the genome (Hu et al. 2015; Bodelon et al. 2016; The Cancer Genome Atlas Research Network 2017; Gagliardi et al. 2020; Warburton et al. 2021; Symer et al. 2022; Fan et al. 2023), most of which contained events that produced fusion transcripts (Fig. 1F). The most frequently integrated loci in our data set were near KLF5/KLF12 (13q22; Supplemental Fig. S3A), MYC (8q24; Supplemental Fig. S3B), and TP63 (3q28; Supplemental Fig. S3C). All of these loci contained four integration events within 1 Mb of the respective gene and all integration events produced HPV–human fusion transcripts. All events in the 8q24 locus were multi-breakpoint events, whereas the 13q22 and 3q28 loci exclusively contained dup-like events (Supplemental Table S2). Another recurrently integrated locus was detected near the repeat-containing ribosomal gene RNA28SN2 (21p11) in four samples (Fig. 1F). Three of the four events were not expressed, but we suspect this was due to a limited ability to map short RNA-seq reads in this low-complexity locus. The one expressed event was a multi-chromosomal event with a fusion transcript detected in a more mappable region overlapping FAM83B on Chromosome 6p12. We compared the expression of MYC, TP63, and KLF5 in samples with and without HPV integration within 1 Mb of the respective genes (Fig. 1G; Supplemental Table S5) and observed that MYC and TP63 expression was significantly increased in the integrated samples, whereas KLF5 expression was not affected (Wilcoxon; MYC P = 0.026; TP63 P = 0.024; KLF5 P = 0.61; Fig. 1G).

HPV structure and copy number heterogeneity within integration events

We defined an HPV integrant as an uninterrupted segment of the integrated virus genome situated between two HPV–human breakpoints. Examples of possible integrant structures are shown in Figure 2A. An HPV breakpoint can be unambiguously linked to another HPV breakpoint when the entirety of the integrant is contained on a single read, and we refer to each unique pairing as a breakpoint pair (Supplemental Table S6). Using such reads, we found evidence for different configurations of viral genome segments—some of which were rearranged and/or concatemerized—between a single breakpoint pair (Fig. 2B). We refer to these as heterologous integrants, and hypothesize that they may arise as a result of unequal amplification of the HPV genome during the replication of HPV–human concatemers, akin to “heterocateny” described by Akagi et al. (2023) and the “onion-skin” structures described by Kadaja et al. (2009) (Fig. 2B). Heterologous integrants were found in 21% (45/212) of breakpoint pairs in our cohort (Fig. 2C; Supplemental Table S6), and up to 15 unique integrant configurations were observed within a single breakpoint pair (Fig. 2C). Figure 2D summarizes the different configurations of HPV integrants present in an example multi-breakpoint event. This sample harbored four connected breakpoints within the HPV genome, splitting it into four segments, which were variably rearranged across the integration event (Fig. 2D). We used read lengths and alignments to distinguish the HPV integrant sizes and structures (Methods), and detected various combinations—some of them periodic—of the four HPV segments between the different HPV–human breakpoint pairs within the event (Fig. 2D).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 2.

Heterologous structures of the HPV genome before and after integration. (A) Schematic of possible HPV integrant structures. The spirals in structures A–E show the portion(s) of the HPV genome that could be contained within an integrant. (B) Schematic showing how several heterologous integrants can exist between a single breakpoint pair, with the size of the integrant varying by n HPV copies. The colors correspond to the regions of the HPV genome as depicted in A. (C) The number of integrant structures between all identified breakpoint pairs within the cohort. Integrants with 2+ identified structures were classified as heterologous. (D) The sizes of the HPV integrants in a multi-breakpoint event with schematics depicting the various integrant structures detected. The HPV genome, in this case, was broken into two segments A and B, and the A segment was further broken into three segments (A.1, A.2, and A.3). These segments were variably rearranged into new structures across the breakpoint pairs. Each point represents the size of an HPV integrant contained on an individual read, which is then grouped together by color (e.g., blue, red, green, light blue) if they do not differ in size by more than 300 bp. Each color in a breakpoint pair thus represents one unique integrant structure, as indicated in the accompanying schematics. (E) The lengths of HPV-aligned reads in four predominantly episomal samples. The existence of HPV episomes and episome concatemers is supported by the accumulation of read counts in bin sizes corresponding to one or more HPV genome copies, as indicated by the dotted lines. (F) Frequency of heterologous integrants in the different integration categories. Only categories harboring heterologous integrants are shown. (G) The percentage of integrants from different HPV types that form single or heterologous structures. (H) The maximum size of the integrant structure in each breakpoint pair in HPV16 and HPV18 integrants, represented as the number of HPV genome copies. (I) Distribution showing the maximum number of HPV copies found in the longest spanning read for each incomplete integrant. The x-axis shows the length of the longest spanning read. Box plots represent the median and upper and lower quartiles of the distribution; whiskers represent the limits of the distribution (1.5 IQR below Q1 or 1.5 IQR above Q3). P-values were calculated using the Wilcoxon rank-sum test and the Fisher's exact test, as indicated in the figure.

Heterogeneity in HPV genome structure has been reported in the episome prior to integration (Rossi et al. 2023). We therefore investigated HPV multimers, i.e., concatenated HPV genomes, and other HPV genome structural variants (SVs) in our four predominantly episomal samples (including one sample with subclonal integration, VAF = 0.00025). HPV genome multimers were found to be the dominant configuration in two samples, while the other two contained predominantly single-copy episomes, although we cannot exclude the possibility that shorter reads may bias these two samples toward detection of smaller configurations (Fig. 2E). Sample HTMCP-03-06-02156 showed an accumulation of read lengths between one and two HPV copies, indicating episomal structural variation, and we indeed found 23 viral SVs in this sample (Supplemental Fig. S4). We thus hypothesize that integrant rearrangement and copy number heterogeneity may partially originate within the episome prior to HPV integration.

Most heterologous integrants (91%) originated from multi-breakpoint events and were not detected in del-like or translocation integration events, compatible with the notion that amplified copy number is associated with their formation (Fig. 2F; Supplemental Table S6). HPV18 integrants were more frequently heterologous than other HPV types (35% vs. 21%; Fisher's exact test, P = 0.0071; Fig. 2G) and contained significantly more viral genome copies per integrant than HPV16 (Wilcoxon, P = 0.040; Fig. 2H).

We were unable to determine the complete integrant for 123 breakpoints because no reads linked the detected breakpoint to any other. For each of these incomplete integrants, we used the longest HPV sequence contained on a read to determine the minimum size of the integrant (Fig. 2I; Supplemental Table S6). Incomplete integrants often contained at least two HPV copies (mean = 2.1) and the longest incomplete HPV integrant was supported by a read that spanned 48 kb (Fig. 2I). In contrast, the largest complete integrant between a breakpoint pair was 37 kb. Thus, even with long-read sequencing, we did not achieve complete resolution of HPV integrant lengths in 123/335 (37%) of HPV breakpoints.

Comparison of two-breakpoint events

Two-breakpoint integration events involving one chromosome comprised 66/129 (51%) of events in our cohort and occurred when one HPV integrant was inserted between two breaks in the human genome. These were associated with either del-like insertions or dup-like expansions (Fig. 1B), which were differentiated based on the alignment patterns of the HPV-containing chimeric reads. Most predicted ecDNAs (7/8) were also two-breakpoint events (Supplemental Table S4). Del-like events were characterized by a copy number loss between the integration breakpoints, in some cases accompanied by amplification of the regions flanking the integration (Fig. 3A). In dup-like events, including six potential ecDNAs, the region between the breakpoints was amplified; in predicted ecDNAs, read alignments were fully contained between and did not extend beyond the breakpoints (Fig. 3A).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 3.

The characteristics of two-breakpoint events and potential ecDNAs. (A) Examples of read coverage patterns from a del-like event, a dup-like event, and a potential ecDNA dup-like event. Orange lines indicate the integration breakpoints. (B) Circular assemblies from three potential ecDNA integration events. The orange portions show the integrated HPV segment, including the viral genes. The direction of HPV gene transcription is shown by a black arrow. The right-most example depicts a complex event in which three nonadjacent human segments have been combined in the potential ecDNA. (C) The size distribution of potential HPV–human hybrid ecDNAs (n = 8). (D) The genomic distance between breakpoints in del-like versus dup-like events. Box plots represent the median and upper and lower quartiles of the distribution; whiskers represent the limits of the distribution (1.5 IQR below Q1 or 1.5 IQR above Q3). (E) The percentage of events occurring in genic (>90% within a gene), intergenic (<10% within a gene), and partially genic (10%–90% within a gene) regions, plotted by integration category. The P-value in D was calculated using a Wilcoxon ranked-sum test. The P-value in E was calculated using a Fisher's exact test.

We predict that candidate ecDNA integration events exist outside the chromosome in a self-contained circular structure, where the human region between the integration breakpoints is bridged by HPV (examples shown in Fig. 3B). The first example shows a near-complete copy of HPV18 connecting two breakpoints from an intergenic region on Chromosome 2 (Fig 3B, left). The second example shows a candidate ecDNA that contains a portion of TP63 along with a portion of HPV45 (Fig. 3B, middle). The third example shows an atypical ecDNA integration event, in which the ecDNA contains part of HPV52 along with a rearranged intergenic segment of Chromosome 13 (Fig. 3B, right). The eight candidate ecDNAs had an average assembly size of 48.8 kb (range 22.4–78.0 kb; Fig 3C) and predominantly involved clade A7 HPV types (4 HPV18, 3 HPV45; Supplemental Table S4). We predicted no HPV16 ecDNAs.

The human–HPV breakpoints in dup-like and del-like integration events were distributed differently across the human genome. Firstly, there was greater spacing between the two breakpoints in dup-like events, with a median distance of 59.1 kb compared to 938 bp (Wilcoxon, P = 7.6 × 10⁻⁶; Fig. 3D; Supplemental Table S7). The dup-like events were also less likely to be completely contained within genic regions and more likely to have a breakpoint within the 3′ or 5′ end of a gene (Fisher's exact, P = 0.040; Fig. 3E; Supplemental Table S7).

HPV integration is associated with full chromosome arm translocations

Translocation integrations, in which two breakpoints occurred on different chromosomes (Akagi et al. 2023), were rare in our study. There were only three instances, involving three different HPV types (HPV18, 52, and 31). The read alignments from the HPV52 translocation integration event are shown in Supplemental Figure S5. This event involved a focal amplification around the site of HPV integration and a segment of HPV containing E6 and E7 (Supplemental Fig. S5A). The Chromosome 8 breakpoint, although within a gene, occurred near the pericentromeric region and was adjacent to a one copy loss upstream of the breakpoint. The Chromosome 1 breakpoint occurred in a genic region on the p-arm and showed copy gain across a ∼70 Mb segment upstream of the integration breakpoint (Supplemental Fig. S5B). This indicates that the Chromosome 1 region was duplicated and recombined with the q-arm of Chromosome 8, the p-arm of Chromosome 8 was lost, and the HPV52 segment was focally amplified around the translocation junction (Supplemental Fig. S5C).

Multi-breakpoint HPV integration events are structurally complex

Multi-breakpoint events, comprising 41/129 (32%) of all events, were the most common type of event in our study. Figure 4A summarizes the number of HPV–human breakpoints per multi-breakpoint event, which ranged from 3 to 32 (Supplemental Table S8). On average, transcribed integration events contained more breakpoints per event than nontranscribed events (Wilcoxon, P = 0.0084; Fig. 4B). Integrations with higher numbers of HPV breakpoints were also associated with more human-only SVs overlapping the integration event (Spearman's correlation, R = 0.69, P = 6.0 × 10⁻⁷) indicating general instability within these regions (Fig. 4C; Supplemental Table S8).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 4.

Complex structural variation is associated with multi-breakpoint integrations. (A) The number of human–HPV and human–human SV breakpoints across multi-breakpoint integration events, and the HPV type and clade in each. (B) The number of breakpoints per event in transcribed and nontranscribed multi-breakpoint events. Box plots represent the median and upper and lower quartiles of the distribution; whiskers represent the limits of the distribution (1.5 IQR below Q1 or 1.5 IQR above Q3). (C) Spearman's correlation between the number of HPV–human breakpoints and the number of human–human SV breakpoints in multi-breakpoint events. (D) Examples illustrating the connectivity between HPV breakpoint pairs in five multi-breakpoint events: three from the MYC-associated locus on Chromosome 8 and the two most rearranged multi-chromosomal events. Dots denote HPV breakpoints along the event, and dotted lines represent the HPV integrants that connect the breakpoints. The dots are colored according to the number of connections that converge at that position in the event. The integrants are colored according to whether single (orange) or heterologous (blue) integrant structures connect the breakpoints. (E) An example breaking down the copy number changes and the proposed structure of an event overlapping MYC. Each color represents a genomic segment on Chromosome 8 in between two human–HPV breakpoints. P-values in B and C were calculated using the Wilcoxon rank-sum test and Spearman's correlation test.

We explored how HPV breakpoints were connected within multi-breakpoint events. Five examples are shown in Figure 4D, three of which originate from the recurrently integrated region on Chromosome 8 near MYC (Fig. 1F). The presence of heterologous integrants between breakpoint pairs is indicated, as is the degree of connectivity of each breakpoint (Fig. 4D). HTMCP-03-06-02058 and HTMCP-03-06-02267 both feature one highly connected breakpoint that pairs with all other breakpoints in the event. This contrasts with the HTMCP-03-06-02149 event, in which each breakpoint connects to two others. The two multi-chromosome events shown in Figure 4D were the two most rearranged events. The HTMCP-03-06-02109 event involved a heavily rearranged Chromosome 3, while Chromosomes 7 and 17 contained one and two breakpoints, respectively. The single breakpoint on Chromosome 7 is connected to almost every other breakpoint on the other two chromosomes. In contrast, HTMCP-03-06-02175 contained three highly connected breakpoints but none connected with all three other chromosomes.

To resolve the copy number patterns that can occur at multi-breakpoint events, we delved deeper into the structure of HTMCP-03-06-02149 event one (Fig. 4E). This event was of particular interest because it encompassed, but did not disrupt, the oncogene MYC. Read depth conversions were used to predict absolute copy number ratios for each segment bridged by an HPV breakpoint pair in the event. The three segments that were connected by nonheterologous integrants, including the region containing MYC, each showed a one copy increase above the baseline adjacent to the event (Fig. 4E). The segment between the heterologous integrant breakpoint pair showed a five-copy increase and this breakpoint pair was associated with heterologous integration, compatible with the notion that amplification may be involved in the generation of heterologous integrants (Fig. 4E). Thus, each of the five copies in the concatenated HPV–human amplification has altered HPV integrant structures. The resulting proposed structure, therefore, contained seven HPV integrants linking the human genome segments, including two uninterrupted copies of MYC, as illustrated in Figure 4E.

HPV integration events show consistent orientation-dependent regulatory patterns

ONT sequencing simultaneously yields DNA methylation information along with DNA sequence data (Simpson et al. 2017). We therefore investigated methylation patterns at HPV integration events, including within the integrated HPV genome and the flanking human regions, and related these to the transcriptional status of the event. Methylation of CpGs within HPV integrants and up to 5 kb on either side of the breakpoints is shown in Figure 5A and B for dup-like and del-like integration events (Supplemental Tables S9, S10). The integration events are oriented according to the direction of HPV transcription. In transcribed events, we observed that human regions upstream of HPV transcription were methylated, while human regions downstream from HPV transcription were unmethylated, and that the HPV late genes (L1 and L2) were methylated, while the early genes (specifically E6 and E7) were unmethylated after the long control region (LCR) up to the 3′ breakpoint (Fig. 5A,B). This pattern was especially pronounced in dup-like events (Fig. 5A). In contrast, events without evidence of HPV expression tended to lack consistent methylation patterns adjacent to HPV and contained partially methylated CpGs across the HPV genome (Fig. 5A,B). In all integrants, the LCR, a noncoding regulatory region that activates HPV early gene transcription (Maki et al. 1996), was invariably hypomethylated compared to the genic region of HPV (Fig. 5A,B).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 5.

The production of HPV fusion transcripts correlates with distinct methylation patterns adjacent to and within the HPV integrant. (A,B) The proportion of HPV-containing reads showing methylation at positions within and adjacent to HPV integrants in (A) dup-like events (amplified region >10 kb) and (B) del-like events. The regions 5 kb upstream of and downstream from (relative to the direction of HPV transcription) are divided into 500 bp bins, and the average methylation frequency for all the CpGs within each bin is shown. Within HPV, the methylation of each CpG bin is shown as a colored dot. The transcriptional status is also indicated for each event (row). All the events are aligned to the start of the genic region (E6 start) for each respective HPV type. The gene model for HPV16 is shown above for general reference. (C) The assemblies of three representative HPV integration events, including their human and HPV gene positions, RNA-seq coverage, and HPV–human fusion junctions. (D) The position of HPV–human RNA fusion points relative to the nearest DNA HPV breakpoint, oriented by the strand of the HPV integrant. The most abundantly expressed junctions (n = 36) are contrasted with all other identified junctions (n = 163). (E) The normalized expression (RPKM) within the 5 kb region downstream from the HPV integrant, stratified by the downstream methylation status. (F) The average downstream methylation (0 to 2500 bp), stratified by HPV transcription status. (G) The normalized expression in reads per kilobase per million (RPKM) within the 5 kb region downstream from the HPV integrant, stratified by the HPV transcription status. (H) Pearson's correlation between downstream methylation and expression in transcribed and nontranscribed events. (I) The difference in the expression (RPKM) upstream of and downstream from (±5 kb) the HPV integrant in events stratified by the downstream methylation status. In all cases, box plots represent the median and upper and lower quartiles of the distribution; whiskers represent the limits of the distribution (1.5 IQR below Q1 or 1.5 IQR above Q3). P-values were calculated using the Wilcoxon rank-sum test and Pearson's correlation, as indicated in the figures. The data points in E–I and the accompanying statistics are from the dup-like and del-like events shown in A and B.

We next investigated how the production of HPV fusion transcripts correlated with methylation patterns across HPV integration events. Integration often truncates or disrupts the HPV genome in such a way that the transcription termination signal for the early genes is lost (normally positioned after E5). We would thus expect transcripts produced from such integrants to read through downstream into the human genome, where transcription might terminate. We show the epigenetic and transcriptional landscapes of HPV integrants using the assemblies of three representative events (Fig. 5C). One event occurred in an intron of NEGR1 upstream of its 3′ exon, and produced an E7^NEGR1 fusion transcript containing the poly(A) tail of the human gene (Fig. 5C). However, the other two events did not overlap any genic elements on the same strand. We observed that, across all samples, HPV–human RNA fusions consistently occurred downstream from the nearest HPV integration breakpoint, especially when focusing on the most abundantly expressed fusion sites (Fig. 5D), and they occurred within the demethylated region adjacent to the HPV integrant (Fig 5C). In fact, events with a demethylated downstream region showed significantly higher normalized expression in the 5 kb bin downstream from HPV compared to events that were methylated (Wilcoxon, P = 2.3 × 10⁻⁵, Fig. 5E). Events that generated HPV fusion transcripts had lower downstream methylation (Wilcoxon, P = 1.9 × 10⁻⁵, Fig. 5F) and higher downstream expression (Wilcoxon, P = 7.3 × 10⁻⁹, Fig. 5G) than nontranscribed events, and expression and methylation in the downstream region were negatively correlated (Pearson's correlation, R = −0.42, P = 0.00086, Fig. 5H). However, the increased expression adjacent to the HPV integrant was unique to the downstream region: in events with a demethylated downstream region, there was a significant imbalance in expression between the upstream and downstream regions, relative to HPV, whereas no such difference was observed in the methylated events (Wilcoxon, P = 1.6 × 10⁻⁵, Fig. 5I).

The methylation patterns in the integrated HPV genomes contrasted with the four samples harboring episomal HPV (Supplemental Fig. S6A), in which viral sequences were hypomethylated compared to the integrated HPV genomes, particularly in the late genic region. The hypermethylation of select HPV genes, therefore, is a common feature of integrated HPV genomes in our cohort and was not seen in samples containing HPV episomes. In the three translocation events, HPV methylation did not follow the standard pattern and we suspect they instead are more influenced by the surrounding human methylation (Supplemental Fig. S6B).

Allele-specific expression and methylation patterns in regions of HPV integration

HPV integrates into only one allele, leaving the other unaffected (Karimzadeh et al. 2023). We leveraged this to explore how HPV integration may have allele-specific impacts on the epigenome and human gene targets on a broader scale. We phased reads into haplotypes and determined the allele from which the HPV-containing reads originated, and then identified differentially methylated regions (DMRs) between the two haplotypes across the tumor genome (Methods). DMR density was compared to a null distribution across the human chromosomes to identify regions with significantly dense clusters of DMRs, referred to as DMR hotspots (Methods). HPV integration sites were then intersected with the DMR hotspots to identify events that potentially affected DNA methylation in an allele-specific manner. Across the samples in our data set, 33 integration event loci overlapped autosomal DMR hotspots (Fig. 6A; Supplemental Fig. S7; Supplemental Table S11). HPV-overlapping DMR hotspots had a median size of 6.1 megabases (Mb; range 1.8–26 Mb) compared to the overall genome-wide median of 3.9 Mb (range 18 kb–56 Mb) and mostly overlapped multi-breakpoint events (Supplemental Fig. S8A,B).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 6.

HPV integration is associated with dysregulation of the methylome and nearby genes on the integrated haplotype. (A) The distribution of DMRs across 10 Mb on either side of HPV integrations in the 33 events that overlapped a DMR hotspot. Each tick represents a DMR. Some event loci were situated <10 Mb from the end of the chromosome or from an unmappable region, resulting in a gap in DMRs before the end of the window (#7, 8, 12, 15, 19, 20, 23, 25, 26, and 29). (B) The direction of methylation changes in the HPV-containing haplotype with respect to the unintegrated haplotype within the phase block containing HPV integration. Adjacent phase blocks are shown as flanking gray bars. Events are ordered identically in A and B. The color of event numbers (center column) indicates the transcriptional status of the event. Events within the red (bottom) and blue (top) boxes show significant DMR enrichment (P adj < 0.05) either unidirectionally (blue box) or bidirectionally (red box) relative to the HPV integration event, as determined by a permutation test of the 500 kb bins flanking the event. (C) The significance of the association between HPV integration and high DMR density at all 147 HPV-integrated regions, using window sizes of 100,000 bp, 500,000 bp, 1,000,000 bp, 2,000,000 bp, and 5,000,000 bp around HPV. (D) The fold change and allele-specific expression (ASE) status of outlier genes (1.5 IQR below Q1 or 1.5 IQR above Q3) within 1 Mb of integration events. The log₂ fold change (log₂FC) of the integrated sample is relative to the median of the cohort. (E) The position of genes with outlier expression (1.5 IQR below Q1 or 1.5 IQR above Q3) relative to sites of HPV integration. Color indicates expression fold change in the integrated sample relative to the median of the cohort. The ASE status, integration event type, and transcriptional status of the event are also indicated. (F) The difference in gene expression fold change (integrated sample/median) at transcribed (yes) and nontranscribed (no) integration events. All box plots represent the median and upper and lower quartiles of the distribution; whiskers represent the limits of the distribution (1.5 IQR below Q1 or 1.5 IQR above Q3). Adjusted P-values in F were calculated using Benjamini–Hochberg-corrected Wilcoxon rank-sum tests. (G) Integrative Genomics Viewer snapshots showing wide (top) and zoomed (bottom) views of the haplotype-specific methylation changes around NR4A3 and HPV integration in HTMCP-03-06-02428 (#31), with reads separated into the two haplotypes (HP1 and HP2). The sample's DMRs, phase blocks, and HPV integration breakpoints are also indicated in the top three tracks. Reads are colored by CpG methylation status, with red indicating methylated and blue indicating unmethylated.

There was variation in DMR distribution at the integration-associated hotspots, in regards to their position relative to HPV and in size and density (Fig. 6A). We therefore tested, in each integration region, if high DMR density was unique to the sample harboring HPV integration by comparing the DMR density in the adjacent 500 kb region to 1000 randomly selected regions (Methods; Supplemental Fig. S9). Three and six events were significantly enriched for DMRs on one (i.e., unidirectional) and both sides of HPV (i.e., bidirectional), respectively, in the integrated sample (P adj < 0.05; Benjamini–Hochberg corrected permutation test; Fig. 6A; Supplemental Figs. S10, S11). Four of the six bidirectional DMR-enriched events produced fusion transcripts (Fig. 6A; Supplemental Table S11), whereas none of the unidirectional DMR-enriched events had HPV–human transcripts detected (Fig. 6A; Supplemental Table S11). The density of DMRs appeared focal around unidirectional DMR-enriched integration events, affecting <1 Mb adjacent to the HPV integrant, while bidirectional enrichment was more distributed, spanning 1–8 Mb from the HPV integrant (Fig. 6A; Supplemental Fig. S12). DMR densities in the remaining 24 events were not statistically different from the same regions in other cervical cancer samples, suggesting the high density of DMRs in these regions may occur independently of HPV integration (Fig. 6A).

We also investigated the direction of DNA methylation changes on the integrated haplotype. All three unidirectional DMR-enriched integration events (events 1–3) were hypomethylated on the HPV-integrated allele, while bidirectional DMR-enriched events showed both hypermethylation (events 31–33) and hypomethylation (events 28–29) (Fig. 6B; Supplemental Table S11). We note that, though methylation changes generally occurred in a uniform direction over a broad region of the human genome in the significantly DMR-enriched events, transcribed events still followed the pattern described in Figure 5 in the CpGs adjacent to the HPV–human breakpoints, suggesting that the local and distal methylation changes may be independent phenomena (Supplemental Fig. S11B).

We also analyzed the association between HPV integration and DMR density across our cohort and observed that regions surrounding HPV integration exhibited a significantly higher density of DMRs compared to 10,000 randomly selected control regions with comparable GC-content (Supplemental Fig. S13A; Methods). This effect was more pronounced closer to the integration event, but was significant up to a 2 Mb (±1 Mb) window (Wilcoxon, P adj < 0.05, Fig. 6C; Supplemental Fig. S13B), indicating that HPV's impact on DNA methylation was strongest near the integration event and gradually diminished around ±1 Mb away from the event (Fig. 6C).

We next searched for nearby genes whose transcription was potentially affected by HPV integration. We limited our analysis to loci within 1 Mb of HPV integration sites, which is approximately the window within which DMR enrichment was observed. We identified 2066 genes near HPV integrants across 68 samples. Of these, 101 genes had outlier expression levels in 31 integrated samples and 43 integration event loci (Fig. 6D,E; Supplemental Table S12; Supplemental Fig. S14). We examined these genes for allele-specific expression (ASE), which would be consistent with a cis-acting element (e.g., HPV integration) altering gene expression on one allele. Of the expression outliers that we could test for ASE, 69% (44/64) showed ASE (Methods; Supplemental Table S12). Fifty percent of the outlier genes (50/101) were within 200 kb of the HPV integration site (in 26 samples and 31 events), and these tended to be overexpressed (48/50 genes; 96%) and have ASE (32/38 tested genes; 84%). Outlier genes between 200 kb and 1 Mb away (51/101) contained significantly fewer that were overexpressed (30/51 genes; 59%; Fisher's exact test, P = 7.4 × 10⁻⁶), indicating that HPV integration may have greater influence on gene expression within 200 kb of the event.

When comparing gene expression in the integrated sample versus the rest of the cohort, genes overlapping expressed HPV integration events (i.e., HPV–human fusion transcripts were detected) showed significantly higher upregulation than those overlapping nonexpressed HPV integration events (Fig. 6F). We also observed greater upregulation of nonoverlapping genes within 200 kb of expressed HPV integration events (Fig. 6F). Genes 200–500 kb away from HPV integration showed no significant difference between expressed and nonexpressed events (Fig. 6F).

The nuclear receptor genes NR4A1 and NR4A3 (Wenzl et al. 2015) were the genes most highly overexpressed in an allele-specific manner from HPV-containing haplotypes (Supplemental Figs. 14, 15A,B). NR4A1 (Chromosome 12) and NR4A3 (Chromosome 9) were detected from two independent integration events in two different samples. Nonintegrated samples showed ASE of NR4A1 and NR4A3 (11 and 2, respectively), but in both cases, the major allele frequency observed in the integrated sample was much higher than the majority of samples (Supplemental Fig. 15A,B). Two outlier genes, NR4A3 and CBARP, overlapped significant DMR-enriched events (event 31 and 29; Fig. 6A,B; Supplemental Table S12), and the remaining seven significant DMR-enriched events were not associated with changes in gene expression within 200 kb, suggesting that allelic methylation and gene expression changes may be two independent consequences of HPV integration.

To determine how HPV may contribute to the upregulation of oncogenes, we explored the DNA methylation landscape around HPV integration in the NR4A3-activated sample. This event was a multi-breakpoint event, with one breakpoint upstream of NR4A3 on Chromosome 9 that connected to a Chromosome 14 locus harboring 12 breakpoints (Supplemental Fig. S16). The Chromosome 9 locus overlapped a DMR hotspot that was hypermethylated on the integrated allele (Fig. 6G, top). However, several DMRs near the NR4A3 promoter were unmethylated on the integrated allele, including an 821 bp DMR in exon two of NR4A3 at the 3′ end of an unmethylated promoter element (Fig. 6G, bottom), distinguishing it from other cervical cancer samples we profiled (Supplemental Fig. S17, bottom). The broad methylation on the HPV-containing allele was also unique to the integrated sample (Supplemental Fig. S17, top). We thus observed, in the NR4A3-activated sample, both HPV-associated large-scale allelic methylation as well as focal demethylation of potential NR4A3 regulatory elements, and high expression of NR4A3 transcripts from the integrated allele.