Integrated HPV16 genomes interact in cis with host chromosomes

As it had previously been shown that interactions between an integrated HPV genome and host chromatin on the same chromosome could lead to changes to host gene expression that could be selected for during carcinogenesis[35], we wished to determine whether these cis interactions were occurring at an earlier stage of cervical neoplasia in our HPV16-positive W12 integrant clones prior to selection. To address this, we first developed a HPV16-specific Region Capture Hi-C protocol (S1A Fig) similar to that published previously[39,40], here using biotinylated RNA baits specific for the HPV16 genome that would select chimeric DNA complexes based upon the virus sequence from a Hi-C library and hence allow determination of any genome interactions in our five integrant clones previously epigenetically characterized[21]. The regions of significant interaction between the HPV16 and host chromatin were determined using GOTHiC software (with significance determined by cumulative binomial test and the threshold set at <0.05) and visualised using the Circos tool (Fig 1). With all clones, there was an absence of interaction reads between deleted regions of the HPV16 genome but also largely across the E1 gene due to the necessary length of RNA bait across this region. In W12 clone G2, reads originate from all MboI fragments across the HPV16 genome, indicating interactions with the host and, although the distribution of the reads was fairly uniform, the greatest percentage of reads came from HPV16 gene E7 (Fig 1A). Significant interactions occurred exclusively with chromosome 5, the chromosome of integration, with the majority likely due to cis interactions with bordering host sequences at the breakpoint. However, upon closer inspection of the single chromosome view (Fig 1A inset), there was a divergence between the majority of reads from the virus—indicating the integration site—and a subset of reads that mapped to a separate region of the host, indicating a long-range interaction between the integrated HPV16 genome in G2 and the host. In clone D2 (Fig 1B), the HPV16 genome also interacted with the site of integration at chromosome 5; however, the virus-host reads, predominantly from the L2 gene, at this scale appear to converge on a single point at the host chromosome (Fig 1B inset). For three further W12 integrant clones (clones H, F and A5), interactions again occurred across the HPV16 genome with the host chromosome of integration, with the majority of reads coming from the E2 portion of the virus genome for all (Fig 1C–1E). In clone H, the HPV16 genome integrated into chromosome 4 and resulted in a large deletion within the host genome (~170 kbp), which is illustrated by the separation of the virus-host reads in the chromosome view (Fig 1C inset). Interestingly, and contrary to our previous publications[15,21,24], we found that W12 clones F and A5 had the same integration site, with virus-host reads converging on the same region of chromosome 4 (Fig 1D and 1E).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. HPV16-specific Region Capture Hi-C determines definitive HPV16 integration sites.

CIRCOS plots show sequence interactions between HPV16 (orange) and host chromosomes (various) for clones (A) G2, (B) D2, (C) H, (D) F and (E) A5. Each line within the circle represents a significant virus-host read indicating an above background interaction between a region of the HPV16 genome and the host. Reads are coloured to match individual HPV16 genes: E6 = green, E7 = orange, E1 = yellow, E2 = blue, E4 = red, E5 = pink, L1 = dark green, L2 = light blue and non-coding regions = black. Percentage of reads coming from different regions of the virus is indicated by the histogram on the outside of the HPV16 genome, which is split into 500 bp windows (red bars). HPV16 RNA bait fragments used in the Capture Hi-C experiment are indicated on the outside of the CIRCOS plot (blue curved lines); interacting reads are largely absent due to either deletion of the genome region during integration or specifically from the E1 gene due to the necessary length of the bait covering this region (see Materials and Methods section for further information). Grey shaded boxes denote the region of HPV16 genome deleted in individual cell lines. Presented data were generated using the Gothic program and plots are not to scale. Insets show zoomed sites of integration, with interaction divergence in clones G2 and H.

https://doi.org/10.1371/journal.ppat.1009875.g001

HPV16 integration site virus-host breakpoint identification at nucleotide resolution

Having developed the probability that virus-host genome interactions did indeed occur soon after HPV16 integration, we next sought to precisely identify the host sequence of the virus-host junctions. To do this with sufficient depth from genomic DNA samples, we developed a novel protocol by which DNA from each integrant clone was enriched for HPV16 sequences along with its flanking regions before sequencing, henceforth known as HPV Integration Site Capture (HISC) (S1B Fig). The resulting sequencing data was analysed for reads mapping to the HPV16 genome with the corresponding human tag being determined. These data were then aligned to the HPV16 genome (S2 Fig) and the human genome (S3–S6 Figs). From all integrant cell lines, sequencing reads mapped to both the HPV16 and host genomes with peaks at only two distinct sites each (S2–S6 Figs, HISC track), regardless of HPV16 genome copy number, demonstrating that only a single 5’ and 3’ breakpoint existed in each W12 clone examined. The finding of only a single 5’ and 3’ for each W12 integrant cloned line assessed was confirmed using ‘LUMPY’, a probabilistic framework for structural variant discovery[41], whereby only four pairs of unique junctions were found with a range of a minimum of 805 split-reads (SRs) (clone H) to a maximum of 6446 SRs (D2). Moreover, the separation of the breakpoint peaks from the HPV16 genome was consistent with termination of RNA-seq reads from separate transcriptome analysis (S2 Fig, RNA-seq tracks) and the known deletion of a proportion of the virus genome in each cell line[21], with the HPV16 transcription profiles additionally being consistent with our previously published quantitative PCR data[15]. Alignment of the breakpoint peaks with host sequence gave a separation distance ranging from ~25–170 kbp across the integrant clones; the separation of peaks here however is consistent with two processes of HPV integration. The greatest distance of 170 kbp seen in clone H occurred due to the deletion of a proportion of the host genome upon ‘direct’ integration of the single HPV16 genome here determined by quantitative PCR (qPCR) of sections of host genome spanning the integration site (S3B Fig). The breakpoint peak separation distance for the other four clones, however, was due to a ‘looping’ mechanism of HPV integration[42], whereby a flanking length of host sequence is amplified during integration of one or more copies of the virus. Again, this was verified through qPCR of host DNA across the integration sites (S4B, S5B and S6B Figs), with transcription clearly occurring through these amplified host regions driven from the integrated HPV16 early promoter (S3A, S4A, S5A and S6A Fig). The resultant genomic effect on the structure of both host alleles is pictured in S3C, S4C, S5C and S6D Figs.

To confirm the virus-host breakpoints determined by HISC, PCR across each junction was carried out and Sanger sequenced, with the verified coordinates of breakage summarized in S7 Fig. Interestingly, a 53 bp region of the E2 gene was found to have been inverted after an intermediate 20 bp deletion, which had not been highlighted by alignment of RNA-seq reads to the HPV16 genome due to bioinformatics processing (S7D Fig). Quantitative PCR (qPCR) was also carried to show that both the 5’ and 3’ breakpoints were unique to each cloned line (other than F and A5), with either very low or non-existent products produced with DNA from an early episomal population of W12 cells (W12par1 p12) or a normal cervix line, NCx/6 (S8 Fig).

Identification of the precise virus-host breakpoints in each cloned line allowed the extent of microhomology at integration sites to be assessed. At all integration sites, at least one end of the insertion involved nucleotides from both genomes directly adjacent to the junction being homologous, with a mode of 5 nt (range 3–5 nt) (S9A–S9D Fig), with clone G2 having 5 nt homology at both ends (S9A Fig). Additionally, when microhomology of 10 nt either side of each breakpoint was compared to that generated from 10,000 random shuffles of each sequence extended to 1,000 nt, 4 of 5 integrant clones (G2, H, F and A5) had statistically significant homology at one flank (S9A, S9C and S9D Fig).

HPV16 integrates into regions of open and transcriptionally active host chromatin

Mapping of the HPV16 integration sites from each cloned cell line allowed investigation of the genomic location and epigenetic landscape into which the virus had inserted (S10 Fig). In 4 of 5 cases, HPV16 had integrated into a gene (D2, TENM2; H, MAPK10; F & A5; RASSF6) with all of these cases arising at introns (S10B–S10D Fig). Despite the integration site in clone G2 occurring intergenically, all insertions occurred at locations of open chromatin, characterized by DNaseI hypersensitivity sites from publicly available normal human epidermal keratinocyte (NHEK) ENCODE datasets (www.encodeproject.org). Additionally, the integration sites showed higher than average levels of histone post-translational modification marks associated with enhancers and transcriptional activity, namely H3K27ac and H3K4me1/2/3, with the marked exclusion of transcriptionally repressed facultative heterochromatin mark, H3K27me3. Noticeably, although integration in clone H did not occur directly into one of these loci, the HPV16 genome is located within 300 kbp of a similar site (S10C Fig). Presence of binding sites of host architectural protein CTCF was found across the sites of integration, although again with higher than average occurrence, and usually within ~20 kbp of the inserted HPV16 genome.

Short- and long-range cis interactions occur between the HPV16 and host genomes

To inspect interactions between the HPV16 and host genomes, Region Capture Hi-C sequencing data was re-visualised using SeqMonk software (Babraham Bioinformatics) and again aligned to NHEK ENCODE datasets at the location of integration for each cloned cell line. Local inspection (~100 kbp window) of interactions between the integrated HPV16 and host genome in clone G2 showed highest levels of interaction corresponding to the 5’ and 3’ virus-host junctions (Fig 2A, red bars), with above background interactions occurring within this region of host genome amplification (Fig 2A, green-orange bars). This profile was consistently found across the Region Capture Hi-C datasets for all clones (Figs 3A and S11), with the exclusion of clone H where the absence of any intermediate interactions between virus:host breakpoints is likely due to deletion of this host region during integration (S11A Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Identification of short- and long-range interactions between integrated HPV16 genomes and the host chromosome in W12 clone G2.

(A) Capture Hi-C data is presented 122.5 kbp across the HPV16 integration locus. The 5’ and 3’ breakpoints of the virus genome (which is not aligned here) are indicated by the tallest red bars and are labelled with black arrowheads, being inverted in comparison to the direction of host sequence due to the ‘looping’ integration mechanism. (B) Capture Hi-C data is presented 700 kbp across the HPV16 integration locus. The black line above the read peaks indicates the genomic window seen in panel A. Peaks of reads indicate regions of the host genome interacting in cis with the integrated virus genome. Short-range interactions between the HPV16 genome and host regions were resolved by consensus using Gothic and are shown beneath the panel, originating from the HPV16 integration site (inverted black arrowhead). (C) Capture Hi-C data is presented 5 Mbp across the HPV16 integration locus. The black line above the read peaks indicates the genomic window seen in panel B. Peaks of reads indicate regions of the host interacting in cis with the integrated virus. Long-range interactions between the HPV16 genome and host regions were resolved by consensus using Gothic and are shown beneath the panel, originating from the HPV16 integration site (inverted black arrowhead). Statistically significant interactions were determined by a cumulative binomial test where adjusted p-value (the q-value), was set a threshold of q<0.05. In each panel, the scale bar represents the normalised read count. Additionally, protein-coding genes are shown in the first track with the direction of each gene indicated by colour (red, forward; blue, reverse), followed by the alignment of ChIP-seq data from the NHEK cell line (ENCODE). Post-translational histone modifications of host enhancers (H3K27ac, H3K4me1; green), active promoters (H3K4me2, H3K4me3; green), repressed chromatin H3K27me3 (red), DNaseI hypersensitivity sites (blue) and CTCF sites (purple) are shown. Coordinates presented for each window are indicated at the top of each figure.

https://doi.org/10.1371/journal.ppat.1009875.g002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Identification of short- and long-range interactions between integrated HPV16 genomes and the host chromosome in W12 clone D2.

(A) Capture Hi-C data is presented 122.5 kbp across the HPV16 integration locus. The 5’ and 3’ breakpoints of the virus genome (which is not aligned here) are indicated by the tallest red bars and are labelled with black arrowheads. (B) Capture Hi-C data is presented 1.4 Mbp across the HPV16 integration locus. The black line above the read peaks indicates the genomic window seen in panel A. Peaks of reads indicate regions of the host interacting in cis with the integrated virus genome. Short-range interactions between the HPV16 genome and host regions were resolved by consensus with Gothic and are shown beneath the panel, originating from the HPV16 integration site (inverted black arrowhead). Statistically significant interactions were determined by a cumulative binomial test where adjusted p-value (the q-value), was set a threshold of q<0.05. In each panel, the scale bar represents the normalised read count. Additionally, protein-coding genes are shown in the first track with the direction of each gene indicated by colour (red, forward; blue, reverse), followed by the alignment of ChIP-seq data from the NHEK cell line (ENCODE). Post-translational histone modifications of host enhancers (H3K27ac, H3K4me1; green), active promoters (H3K4me2, H3K4me3; green), repressed chromatin H3K27me3 (red), DNaseI hypersensitivity sites (blue) and CTCF sites (purple) are shown. Coordinates presented for each window are indicated at the top of each figure.

https://doi.org/10.1371/journal.ppat.1009875.g003

Upon increasing the window range across the integration loci to ~700 kbp (with re-normalised read depth), peaks of genome interaction outside of the initial 100 kbp window could be seen in clone G2 (Fig 2B) and clone D2 (Fig 3B). Multiple short-range (<500 kbp) interactions were present ranging up to ~240 kbp and ~360 kbp away from the site of integration in clones G2 and D2, respectively (Figs 2B and 3B). As before, LUMPY analysis of HISC sequencing libraries confirmed that true chimeric DNA reads were only found at virus integrant breakpoints and, therefore, with a very high likelihood (binomial test, p<0.00001) that virus:host looping interactions were not aberrantly called due to covalent linkage between HPV16 and host genomes at further breakpoints. Interestingly, sites of high intensity of interaction with the host genome in both cell lines (Figs 2B and 3B, red bars) overlapped with histone marks of transcriptional activation or enhancers and DNaseI hypersensitivity sites, whereas these sites only correlated with host CTCF binding sites in clone G2 (Fig 2B, purple dashes). The presence of a long-range interaction (>500 kbp) was seen at this scale in clone D2 ~530 kbp from the integration site with the 3’ end of host gene TENM2 (Fig 3B). However, long-range interactions in clone G2 were only visible when expanding the window range to 5Mb, whereupon two clear interactions were determined downstream of the integration site, with the furthest ~900 kbp from the HPV16 integration site (Fig 2C). This interaction was located at Chr5:53,520,000 within the first intron of host gene ARL15, coinciding with a cluster of host CTCF binding sites at a topologically associating domain (TAD) boundary (Fig 2C, purple dashes) and is also the long-range interaction visible on the CIRCOS plot (Fig 1A, inset). Additionally, to statistically confirm that the aggregation of interactions between virus and host genomes could not have occurred by chance, we used a hypergeometric test (Fisher’s exact test) to compare the region of interactions in both clones G2 and D2 (Chr5:52,000,000–54,000,000 and Chr5:167,000,000–167,600,000, respectively) to equivalent length, adjacent ‘blank’ regions (Chr5:50,000,000–52,000,000 and Chr5:166,400,00–167,000,000 for clones G2 and D2, respectively) (Figs 2C and 3B, respectively) that had a level of background reads analogous to those seen across the chromosome, which clearly showed a significant difference (p<0.00001) in both cases.

To validate that direct interaction might occur between the integrated HPV16 genome and the ARL15 intron in clone G2, three-dimensional (3D) DNA fluorescence in situ hybridisation (FISH) was carried out (Fig 4). Three fluorescent DNA probes were produced to hybridise to either the integrated HPV16 genome (green), the ARL15 site of interaction (red), or a control region of the host genome (purple) with the same linear distance in the opposite direction (Fig 4A). Only cells containing one HPV16 signal and two copies of both the control and ARL15 probes were analysed. A representative image from the resulting dataset (30,000 cells) is shown in Fig 4B and analysis of the 3D distances (x, y and z planes) indicated that, on the integrated chromosome, the HPV16 and ARL15-specific probes were significantly closer together than the HPV16 and control-specific probes (Fig 4C and 4D). Additionally, it was found that the distance between the control and ARL15-specific probes on the integrated allele was significantly lower than that on the unintegrated allele (Fig 4E and 4F). In comparison, the HPV16 signal in W12 clone H was rarely found near to that of either ARL15 or control region, consistent with the HPV16 integration site being on a different chromosome, whilst HPV16 signal was never found in HPV-negative Raji cells (S12 Fig). Therefore, collectively DNA FISH and capture Hi-C data strongly support a direct interaction between the integrated HPV16 genome and the host ARL15 gene in clone G2 (Fig 4B, inset).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Validation of HPV16-host genome cis interactions in W12 clone G2 by fluorescence in situ hybridisation (FISH).

(A) Schematic detailing the complementarity of the DNA probes used on the integrated and unintegrated alleles of a portion of chromosome 5 (51–54 Mbp) in W12 clone G2 to confirm interaction between the HPV16 genome (black arrow) and ARL15 gene (red arrow): Control probe (51,676,020–51,873,551; purple), HPV16 probe (green), and ARL15 probe (53,473,886–53,584,235; red). Possible interactions between probe regions are also highlighted. (B) Representative images of the probes hybridised to W12 clone G2 genome of one cell in a 3D FISH experiment (nucleus boundary, blue), a composite image with DAPI (NUC, blue) stain, and interpretation of the associated chromosome spatial conformations. A zoom of the interaction signals is inset. (C) Box-whiskers plot and (D) frequency distribution chart of the distance between both sets of FISH probes in the integrated allele of chromosome 5: HPV16:ARL15 (red box) and HPV16:control (purple box). (E) Box-whiskers plot and (F) frequency distribution chart of the distance between the Control and ARL15 probes in both the integrated (green) and unintegrated (grey) alleles. Lower and upper whiskers denote the 10th and 90th percentiles, respectively, of the distribution. The lower and upper limits of the boxes indicate the 25th and 75th percentiles, respectively. Solid line in the box denotes the median. Numbers below the box plots denote mean ± SEM (n = 585) from which an unpaired, two-tailed Students T-test was conducted; *p<0.05, ****p<0.0001.

https://doi.org/10.1371/journal.ppat.1009875.g004

HPV16 genome integration sites with short-range (~25kbp) host genome amplification do not affect host topologically associating domain boundaries

Since interactions between integrated HPV16 and host genomes had been confirmed, we next asked whether these interactions could lead to structural changes in local nuclear architecture. Hi-C libraries, containing global host:host interactions, from both clones G2 and D2 were sequenced and, due to the different sites of integration on chromosome 5, these datasets were compared against each other using Insulation Score analysis (Fig 5). This methodology determines how insulated (i.e. isolated from interaction in three-dimensional space) one genomic region is from its neighbouring regions[43] and, as such, results in an albeit lower resolution but more robust analysis of genome interactions across, for example, a 5Mb window. Consequently, regions with the most positive insulation score can be defined as topologically associating domain (TAD) boundaries. Heatmaps of clone G2 Hi-C data (Fig 5A) showed an interaction profile consistent with data from clone D2 (Fig 5B), with interactions occurring across the 5Mb window including clearly defined regions of interacting host DNA up to ~1Mb, approximating to the size of an average mammalian TAD. Insulation score analysis also determined across the 5Mb window no significant change between the two interaction profiles despite the presence of HPV16 genomes at the integration site in clone G2 (Fig 5C). Indeed, this finding was replicated when 2.5Mb either side of the clone D2 HPV16 integration site was compared to that of clone G2 (Fig 5D–5F). Hence, HPV16 genome integration and interactions with the host chromosomes did not appear to be affecting the local nuclear architecture at the hierarchical level of TADs and their boundaries.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. HPV16 genome integration does not disrupt host topologically associating domain (TAD) boundaries at the integration site in W12 clones G2 and D2.

(A, D) Hi-C data for W12 clone G2 is compared to (B, E) Hi-C data for W12 clone D2 using (C, F) an insulation score interpreting the boundaries of topologically associating domains (TADs) for HPV16 integration sites within W12 clone G2 (Chr5: 50–55 Mbp; left column) and W12 clone D2 (Chr5: 164.6–169.6 Mbp, right column) showing no significant change to either window. Black arrowhead = integration site.

https://doi.org/10.1371/journal.ppat.1009875.g005

HPV16 genome integration does modulate intra-TAD host:host genome interactions and host gene expression

As HPV16 genome integration, and resulting interactions with the host genome, did not appear to affect TAD positioning, we sought to determine how far within the host:host chromosome interactions ‘loops’ from the HPV16 genome occurred. To the existing Region Capture Hi-C datasets, we aligned consensus W12 TAD boundaries (Figs 6 and 7), determined from total Hi-C libraries of clones D2 and G2, which showed no difference in TAD calling between clones regardless of HPV16 genome integration, confirming the finding of our insulation score analysis (Fig 5). We found that all HPV16:host genome interactions occurred within, or at the boundary of, the TAD of integration in both clones (Figs 6A and 7A), and all loops in clone D2 occurred exclusively with the TENM2 gene (Fig 7A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. HPV16 genome integration and virus:host genome interactions lead to significant modulation of host gene expression in W12 clone G2.

(A) Capture Hi-C data is presented across Chr5:50–55 Mbp for clone G2. HPV16 integration site is indicated with a black arrowhead and CTCF sites (purple) aligned across the top of the panel. Aligned protein coding genes are shown in the top track (rightward, red; leftward, blue) with the extent of W12 topologically associating domains (TADs) shown below. Charts indicating the transcript level of host protein coding genes within the 5 Mb region of clone G2 relative to (B) W12 episomal level and (C) mean control level across all other integrant clones. All data is shown as a Log2 fold-change with statistically significant changes indicated by green bars (p<0.05, negative binomial Wald test). Gene length is indicated by width of the corresponding bar.

https://doi.org/10.1371/journal.ppat.1009875.g006

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. HPV16 genome integration and virus:host genome interactions lead to significant modulation of host gene expression in W12 clone D2.

(A) Capture Hi-C data is presented across Chr5:164.4–169.6 Mbp for clone D2. HPV16 integration site is indicated with a black arrowhead and CTCF sites (purple) aligned across the top of the panel. Aligned protein coding genes are shown in the top track (rightward, red; leftward, blue) with the extent of W12 topologically associating domains (TADs) shown below. Charts indicating the transcript level of host protein coding genes within the 5 Mb region of clone D2 relative to (B) W12 episomal level and (C) mean control level across all other integrant clones. All data is shown as a Log2 fold-change with statistically significant changes indicated by green bars (p<0.05, negative binomial Wald test). Gene length is indicated by width of the corresponding bar.

https://doi.org/10.1371/journal.ppat.1009875.g007

Next, we addressed whether HPV16 genome integration and interactions with the host chromosome led to any changes in local host gene expression. Transcribed RNAs across a 5Mb window spanning the HPV16 integration site were assessed through comparison of total transcriptome RNA-seq data from an individual integrant clone with transcript levels in both the parental episomal W12 cell line and, separately, an average of all other available datasets from W12 integrant clones with a different integration site. Analysis of protein coding RNAs from clone G2 found that host genes were both up- and down-regulated across the integration locus and in some cases unchanged, with no individual direction of change (i.e. up-regulation or down-regulation) of differential expression restricted to individual TADs (Fig 6B and 6C); a consistent finding across all cloned cell lines analysed. Despite all but one change being less than ±2-fold here (a 2.64-fold decrease in ITGA2 transcript level in comparison to episomal cells (Fig 6B)), some genes were found to be significantly down-regulated versus both episomal and mean integrant clone levels (Fig 6B and 6C, respectively; green labelled genes) (PARP8, ITGA2, DHX29, SKIV2L2) whilst, interestingly, the only significantly up-regulated gene was that with a confirmed HPV16:host interaction, ARL15 (1.23-fold; p<0.05) in comparison to mean integrant clone levels (Fig 6C). Compellingly, upon comparative analysis of Hi-C libraries between clones G2 and D2 (a more direct comparison of the local interactions than Insulation Score analysis), a decrease in a host:host interaction was found within the TAD of integration in clone G2 (S13A Fig, yellow outlined triangle), in which bounds aligned with the HPV16:host interaction at the TAD boundary within the ARL15 gene and just downstream of ITGA2. Hence, some changes of host gene expression within TADs were found to correlate with modulated host:host interactions as well as novel virus:host interactions.

Analysis of host gene expression at the HPV16 integration site of clone D2 again found that transcription was both up- and down-regulated across the 5Mb region, whilst some genes appeared unaffected (Fig 7). The host gene into which HPV16 had integrated and exclusively interacted with via looping, TENM2, was down-regulated in comparison to episomal cells (Fig 7B) likely due to the very high transcript level in parental W12Par1 cells (S14 Fig). However, in comparison to the integrant clone average, TENM2 transcript was up-regulated 4.79-fold (Fig 7C). This highlights the importance of using the mean of several cloned lines as a comparate for modified host transcript levels, in comparison to true low passage W12 episomal cells, as these cloned lines were isolated from the parental population using at least 60 population doublings. No significant changes in host:host interaction were found here within the TAD of HPV16 integration or the adjacent TADs (S15 Fig). Interestingly, all other clones that exhibited HPV16 integration within a host gene showed up-regulation of expression of that gene in comparison both to episomal and integrant clone lines. Expression of RASSF6 in clones F and A5 increased in both cases ~ 2.3-fold versus episomal cells (Fig 8B and 8D, respectively) and ~1.6-fold in comparison to mean integrant levels (Fig 8C and 8E, respectively), with p<0.0001 in all cases. This was in contrast to the varied up- and down-regulation of expression across the 5Mb window when comparing all other analysed genes in clones F and A5 (Fig 8B–8E, respectively).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. HPV16 genome integration leads to significant, but differential, modulation of host gene expression in W12 clones F and A5.

(A) Representative Capture Hi-C data is presented across Chr4:72–77 Mbp for clone F and A5. HPV16 integration site is indicated with black arrowheads (breakpoints) and CTCF sites (purple) are aligned across the top of the panel. Aligned protein coding genes are shown in the top track (rightward, red; leftward, blue) with the extent of W12 topologically associating domains (TADs) shown below. Charts indicating the transcript level of host protein coding genes within the 5 Mb region of clone (B-C) F and (D-E) A5 relative to (B, D) W12 episomal level and (C, E) mean control level across all other integrant clones. All data is shown as a Log2 fold change with statistically significant changes indicated by green bars (p<0.05, negative binomial Wald test). Gene length is indicated by width of the corresponding bar.

https://doi.org/10.1371/journal.ppat.1009875.g008

In clone H, despite deletion of some of the host coding exons, expression of MAPK10 was increased 5.06-fold and 4.47-fold (p<0.001 in both cases) versus W12 episomal and mean integrant transcript levels, respectively (Fig 9B and 9C, respectively). Moreover, further interrogation of RNA-seq data determined chimeric HPV16:host RNA-seq reads indicating both breakpoint fusion transcripts and spliced transcripts from the integrated HPV16 genome (including differential HPV16 exon expression) into an adjacent host exon (ENSE00001811960; Chr4:86,952,584), which was the likely cause of the overall increase in expression levels of this gene (S16 Fig). Indeed, splicing events from the integrated HPV16 genome into the host was seen across three of the five integrant clones analysed (S1 Table), including clone G2 where intergenic HPV16 integration led to spliced fusion transcripts with non-coding host DNA (verified through PCR and Sanger sequencing), presumably through cryptic splice acceptor sites (S17 Fig). Interestingly, all splicing events occurred with host DNA within the region of host DNA amplification flanking the HPV16 integration site in clone G2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. HPV16 genome integration leads to significant modulation of host gene expression in W12 clone H.

(A) Capture Hi-C data is presented across Chr4:84.5–89.5 Mbp for clone H. HPV16 integration site is indicated with a black arrowheads (breakpoints) and CTCF sites (purple) aligned across the top of the panel. Aligned protein coding genes are shown in the top track (rightward, red; leftward, blue) with the extent of W12 topologically associating domains (TADs) shown below. Charts indicating the transcript level of host protein coding genes within the 5 Mb region of clone H relative to (B) W12 episomal level and (C) mean control level across all other integrant clones. All data is shown as a Log2 fold-change with statistically significant changes indicated by green bars (p<0.05, negative binomial Wald test). Gene length is indicated by width of the corresponding bar.

https://doi.org/10.1371/journal.ppat.1009875.g009

HPV16 genome integration modulates host gene expression across the chromosome

To further investigate the effect of HPV16 integration on host gene expression, the variance in gene expression in the genomic regions adjacent to the HPV16 integration site was compared with that of the whole chromosome. Here, comparing RNA-seq data of 100 genes either side of the integration site (coordinates of genome ranges given in S2 Table) within an individual clone to an average of all other available datasets from W12 integrant clones with a different integration site, host transcript level variance (from both protein coding and non-coding regions) was calculated by grouping five adjacent genes into ‘bins’ (average genomic range of 1.1Mb ±0.3Mb) and then comparing the range and variance of gene expression within each bin to that of a mean level from across the whole chromosome (Fig 10). Although individual protein-coding genes at HPV16 integration sites had previously shown statistically significant changes in comparison to the integrant clone average (Figs 6C, 7C, 8C, 8E and 9C), averaging transcript levels from both protein-coding genes and non-coding genes led to statistical significance at only the integration site in clone A5 (Fig 10E and 10J). Expression of host genes across all regions indeed appeared variable, however bin mean fold-changes of greater than or equal to 2 in comparison to chromosome mean expression level were limited to only one bin in all clones (Fig 10A–10C and 10E: G2, bin 9; D2, bin -18; H, bin 5; A5, bin 0) except clone F with two adjacent bins (Fig 10D: F, bins -7 and -6). Interestingly, all bins with mean fold-changes of greater than or equal to 2 consisted of up-regulated mean expression (Fig 10A–10D), apart from clone A5 exclusively with down-regulated mean expression at the site of integration (Fig 10E) despite the increase in RASSF6 transcript levels (Fig 8E). Bins at, and adjacent to, the HPV16 integration site were also found to have highly significant gene expression variance (p<0.01), although this finding was restricted to integrant clones without virus:host cis interactions (Fig 10H–10J). Statistically significant variance was found as far away from the HPV integration site as 5.6Mb in clone H (Fig 10H, bin +5) and 6.4Mb in clone F (Fig 10I, bin -6). Hence, our findings suggest that integration of HPV16 genomes into host chromosomes could have direct modulatory effects on host gene expression beyond both the immediate locus of integration and also the indirect effects of virus oncogene expression on average chromosomal expression.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Variance in host gene expression across the host genomic regions containing the HPV16 integration site in W12 clones.

Each left panel indicates the range and variance of host gene expression in W12 integrant clones (A) G2, (B) D2, (C) H, (D) F and (E) A5, focussing on 100 genes either side of the HPV16 integration site. For each clone, host transcript expression levels were compared with the mean of all other integrant clones. In each panel, the HPV16 integration site is centred on ‘bin 0’. Each bin contains five genes, with no overlap between bins. The box and whisker plots illustrate the range of gene expression levels within each bin, with the black bar indicating median values, the box the interquartile range (IQR) and the whiskers the range. The mean gene expression across the whole chromosome is indicated by the solid blue line, while the mean level of gene expression across individual bins is shown by the dotted blue line. The mean variance of gene expression across the whole chromosome is indicated by the solid red line, while the mean variance within individual bins is shown by a dotted red line. All transcript data is Log2 transformed. Each right hand panel shows the significance of the variance (F-test) in gene expression within each bin for clones (F) G2, (G) D2, (H) H, (I) F and (J) A5. Each point represents a five-gene bin, corresponding to those in the left-hand panels. The horizontal lines indicate the significance of the variance in each bin, compared with the variance in gene expression across the whole chromosome (above the dashed red line, p<0.05; above the dashed pink line, p<0.01).

https://doi.org/10.1371/journal.ppat.1009875.g010

Cell culture

Detailed descriptions of the W12 system have been published previously[10,74,75], including W12 integrant clone generation[15,24] from the parental, episomal W12 culture series-2 (W12Ser2)[76], whereby HPV16 integrant clones used here more likely represent those that occur naturally after HPV16 infection due to non-competitive selection, rather than interferon-induced selection pressure used previously[37]. As such, the five W12 clones used here were episome-free, did not express the HPV16 transcriptional regulator E2[15] and were grown in monolayer culture in order to restrict cell differentiation and maintain the phenotype of the basal epithelial cell layer[77]. Additionally, W12 clones were analysed at the lowest possible passage after cloning (typically p3 to p8) in order to minimise any effects of genomic instability caused by deregulated HPV16 oncogene expression[21]. Furthermore, the initial episomal W12 (Par1) cell line[8] and HPV-negative normal cervix line, NCx/6[10,15], were used as additional controls.

HPV16-host breakpoint and splice junction confirmation

To verify chimeric DNA sequences of HPV16-host breakpoints determined by Capture-seq, and to confirm splice junction sequences from clone G2 RNA-seq analysis, primers were designed using either Primer3 (Primer3Web) alone or Primer-BLAST (NCBI) specific to the DNA sequence (or cDNA sequence from reverse transcribed clone G2 RNA samples (QuantiTect Reverse Transcription Kit, Qiagen)) for PCR (PCR SuperMix High Fidelity, Thermofisher). PCR products were gel extracted and then Sanger sequenced using both 5′- and 3′-end primers to confirm reads from each end of the product. Each analysis was carried out in duplicate. Primer pairs used for PCR of HPV16-host breakpoints and clone G2 splice junctions are given in S3 and S4 Tables, respectively.

qPCR of HPV16-host breakpoints and genomic DNA

Primers were designed using either Primer3 alone or Primer-BLAST (NCBI) specific to the chimeric DNA sequence of HPV16-host junctions determined by Capture-seq to verify integration sites, as well as host DNA spanning integration sites to determine copy number after HPV16 integration. Primers used for qPCR are given in S5 and S6 Tables. Host DNA copy number was quantified by comparison to TLR2 and IFNβ, as reported previously[15]. Conditions used for all primer pairs on an Eppendorf Mastercycler Realplex were: 95°C for 2min; 40 cycles of 95°C for 15sec, 58°C for 20sec, 72°C for 15sec, 76°C for 5sec and read; followed by melting cure analysis from 65°C to 90°C to confirm product specific amplification.

DNA Fluorescence In-Situ Hybridisation (FISH)

BAC clones RP11-467N14 (control locus) and CTD-2015C9 (ARL15 locus) were purchased from (Thermofisher), whereas the HPV plasmid pSP64-HPV16 was prepared in house. BAC/plasmid DNA was purified using the NucleoBond BAC100 kit (Macherey-Nagel), and labelled with aminoallyl-dUTP by nick translation. After purification, 0.5–1 μg labelled BAC DNA was coupled with Alexa Fluor 488, Alexa Fluor 555 or Alexa Fluor 647 reactive dyes (Life Technologies) according to the manufacturer’s instructions, and DNA FISH was performed and analysed as described elsewhere[78,79].

Chromatin crosslinking

Formaldehyde crosslinking of 30 million cells was performed by supplementing standard EGF positive culture medium with formaldehyde to a final concentration of 2% and was carried out for 10 min at room temperature. Crosslinking was quenched by the addition of ice-cold glycine to a final concentration of 125 mM. The adherent cells were scraped from the cell culture plates after crosslinking, collected by centrifugation (400 g for 10 minutes at 4°C), and washed once with PBS (50 ml). After centrifugation (400 g for 10 minutes at 4°C), the supernatant was removed, and the cell pellets were snap-frozen in liquid nitrogen and stored at -80°C.

Hi-C library generation

Cells were thawed on ice, and then lysed on ice for 30 minutes in 50 ml freshly prepared ice-cold lysis buffer (10 mM Tris-HCl pH 8, 10 mM NaCl, 0.2% Igepal CA-630, one protease inhibitor cocktail tablet (Roche complete, EDTA-free)). Following the lysis, nuclei were pelleted (650 g for 5 minutes at 4°C), washed once with 1.25 x NEBuffer 2, and then re-suspended in 1.25 x NEBuffer 2 to make aliquots of 5–6 million cells for digestion. SDS was added (0.3% final concentration) and the nuclei were incubated at 37°C for one hour (950 rpm). Triton X-100 was added to a final concentration of 1.7% and the nuclei were incubated at 37°C for one hour (950 rpm). Restriction digest was performed overnight at 37°C (950 rpm) using 800 units MboI (NEB) per 5 million cells. Restriction fragment ends were filled in using Klenow (NEB) with dCTP, dGTP, dTTP and biotin-14-dATP, and the blunt-ended DNA was ligated following the in-nucleus ligation protocol described previously[80], with minor modifications. Prior to ligation, excess salts and enzymes were removed by centrifugation (600 g for 5 minutes at 4°C) and the cell pellet was re-suspended in 995 μl of 1 x ligation buffer (NEB) supplemented with BSA (100 μg/mL final concentration). The ligation was carried out using 2000 units of T4 DNA ligase (NEB) per 5 million aliquot of cells, at 16°C for 4 hours, followed by 30 min at room temperature. Chromatin was then de-crosslinked overnight at 65°C in the presence of proteinase K (Roche), purified by phenol and phenol-chloroform extractions, precipitated with ethanol and sodium acetate and re-suspended in TLE (10 mM Tris-HCl pH 8.0; 0.1 mM EDTA). The DNA concentration was measured using the Quant-iT PicoGreen assay (Life Technologies). 40 μg of Hi-C library DNA were incubated with T4 DNA polymerase (NEB) for 4 hours at 20°C to remove of biotin from non-ligated fragment ends, followed by phenol/chloroform purification and DNA precipitation overnight at -20°C. DNA was sheared to an average size of 400 bp using the Covaris E220 (settings: duty factor: 10%; peak incident power: 140W; cycles per burst: 200; time: 55 seconds). End-repairing of the sheared DNA (using T4 DNA polymerase, T4 DNA polynucleotide kinase, Klenow (all NEB)) was followed by dATP addition (Klenow exo-, NEB) and a double-sided size selection using AMPure XP beads (Beckman Coulter) to isolate DNA ranging from 250 to 550 bp. Biotin-marked ligation junctions were immobilised using MyOne Streptavidin C1 Dynabeads (Invitrogen) in binding buffer (5 mM Tris-HCl pH 8.0, 0.5 mM EDTA, 1M NaCl) and after stringent washing in the same buffer at 55°C for 10 min ligated to the custom capture-seq adapter using 1600 units of T4 DNA ligase (NEB) for 2 hours at room temperature. The immobilised Hi-C libraries were amplified using the custom primers PE PCR 1.0.33 and PE PCR 2.0.33 with 7–9 cycles. After PCR amplification, the Hi-C libraries were purified with AMPure XP beads (Beckman Coulter). Quantity and integrity of the Hi-C libraries was determined by Bioanalyzer profiles (Agilent Technologies).

Genomic DNA library generation

Cells were thawed, lysed and nuclei were isolated as described above. Nuclei from 5–6 million cells were treated with SDS and Triton X-100 as described for the generation of Hi-C libraries. All Hi-C specific steps, such as MboI digestion, restriction fragment end fill-in, blunt end ligation and the removal of biotin from un-ligated restriction fragment ends were mock performed by replacing the respective enzymes with an equal amount of water. All other steps were performed as described for the generation of the Hi-C libraries. The biotin-streptavidin pull down was omitted and, therefore, the ligation of the custom sequence adapters was done in solution by adding 4 μl adaptors (30 μM) and 1600 units T4 DNA ligase (NEB). The ligation was carried at for 2 hours at room temperature on a rotating wheel in 1x ligation buffer (NEB). Pre-capture PCR amplification was carried out using the custom primers PE PCR 1.0.33 and PE PCR 2.0.33 with 7–8 cycles. The amplified libraries were purified with AMPure XP beads (Beckman Coulter) and their quantity and the quality were assessed by Bioanalyzer profiles (Agilent Technologies).

Capture RNA bait library design

120-mer capture RNA baits were bioinformatically designed to both ends of MboI restriction fragments overlapping the HPV16 genome. Requirements for target sequences were as follows: GC content between 25% and 65%, no more than two consecutive Ns within the target sequences, and maximum distance to a MboI restriction site 330 bp. For short MboI fragments, where 120-mer RNA baits originating from both ends would have overlapped (potentially interfering with optimal hybridization to Hi-C libraries), only the Watson (coding or sense) strand was used for capture RNA bait design, and if necessary the baits were trimmed to minimum length no shorter than 97 nt. This resulted in the design of 16 RNA bait sequences (S7 Table) covering the MboI restriction fragment ends of the entire HPV16 genome, with the exception of two fragments too short (18 and 63 bp, respectively) for capture RNA bait design.

Biotinylated RNA bait library for Region Capture Hi-C generation

The process of HPV16-specific Region Capture Hi-C was carried out essentially as previously published[81]. DNA sequences encoding for the 16 RNA bait sequences, with different restriction enzymes sites on each side (BglII on one site and either HindII or SpeI on the other), which were separated by a 3 bp random spacing sequence, were ordered as two gBlocks (Integrated DNA Technologies) and cloned into plasmid vectors using the Zero Blunt TOPO cloning kit with One Shot TOP10 Chemically competent cells according to manufacturer’s instructions. Both gBlocks were extracted from plasmid DNA by EcoRI (30 units) restriction enzyme digestion at 37°C for 2 hours. Having a BglII and another restriction enzyme sites on the other end, enabled BglII side specific ligation of a T7 promoter sequence adapter with a BamHI overhang essentially as per manufacturer’s instructions, preventing the generation of overlapping complementary transcripts. These adapters were generated by annealing T7_promoter_adapter_1 and T7_promoter_adapter_2 as per manufacturer’s insctructions. Digestion with both restriction enzymes and adaptor ligation were done in one reaction simultaneously, in the presence of BamHI. Two reactions containing 700 ng of gBlock1 DNA or 850 ng of gBlock2 DNA, 30 units BglIl each, 100 units BamHI each, 5-fold molar access of pre-annealed T7 promoter adapters and either 80 units HindIII (NEB) or 40 units SpeI (NEB) were incubated at 37°C for 2 hours in 1x T4 DNA ligase buffer (NEB). Following this incubation 1200 units T4 DNA ligase (NEB) were added to each reaction and incubated at 25°C for 3 hours. The samples were then run on a 1% agarose gel and specific bands at 180 bp were cut out and gel purified. Equimolar amounts of sequences were in vitro transcribed according to manufacturer’s instructions using the T7 MegaScript kit (Ambion) with biotin-labelled UTP (Roche). The RNA was then purified using the MEGAclear kit (Ambion) following the manufacturer’s instructions.

Biotinylated RNA bait library generation for HISC

Four consecutive and non-overlapping fragments from the pSP64 HPV16 bacterial artificial chromosome were PCR amplified using the Expand High Fidelity PCR system (Roche). This resulted in complete coverage of the HPV16 genome. T7 promoter sequences (Roche) were added to one side of the PCR product during this PCR amplification, enabling subsequent directional in vitro transcription. Sequences were in vitro transcribed in the presence of biotin-UTP and purified as described above, followed by fragmentation to 120 nt with 4 nM magnesium chloride at 95°C for 7 min. Fragmented biotinylated RNA was purified by Isopropanol precipitation followed by two subsequent washes with 75% ethanol.

Solution hybridization Region Capture of Hi-C libraries

500 ng to 2000 ng of Hi-C library DNA or genomic DNA library were concentrated using a vacuum concentrator (Savant SPD 2010, Thermo Scientific) and then re-suspended in 5 μl dH₂O. 2.5 μg mouse cot-1 DNA (Invitrogen) and 2.5 μg sheared salmon sperm DNA (Ambion) were added as blocking agents. To prevent concatemer formation during hybridization, 1.5 μl blocking mix (300 μM) was added (equimolar mix of P5_b1_for_33, P5_b1_rev_33, P7_b2_for and P7B2_rev (see list X for sequences)). Biotinylated RNA baits were used in a ratio 1:12 to Hi-C libraries (25 ng biotinylated RNA baits per 300 ng of Hi-C library) and supplemented with 30 units SUPERase-In (Ambion). Biotinylated RNA baits for capture DNA-Seq were used in a ratio of 1:3.33 (300 ng RNA baits per 1,000 ng genomic DNA library) and supplemented with 30 units SUPERase-In. The DNA was denatured at 95°C for 5 min in a PCR machine (PTC-200, MJ Research; PCR strip tubes (Agilent 410022)) and then incubated with the biotin capture RNA at 65°C, in hybridization buffer (5 x SSPE (Gibco), 5 x Denhardt`s solution (Invitrogen), 5 mM EDTA (Gibco), 0.1% SDS (Promega)) for 24 hours, in a total reaction volume of 30 μl. Captured DNA/RNA hybrids were enriched using Dynabeads MyOne Streptavidin T1 beads (Life Technologies) in binding buffer (1 M NaCl, 10 mM Tris-HCl pH 7.5, 1 mM EDTA) for 30 minutes at room temperature. After washing (once in wash buffer 1 (1 x SSC, 0.1% SDS) for 15 minutes at room temperature, followed by three washes in wash buffer 2 (0.1 x SSC, 0.1% SDS) for 10 minutes each at 65°C), the streptavidin beads (with bound captured DNA/RNA) were re-suspended in 30 μl 1 x NEBuffer 2. Post-capture PCR amplification was carried out using between six to nine cycles using primer pairs that consisted of one TruSeq adapter reverse compliment and the TruSeq universal adapter from streptavidin beads in multiple parallel reactions, which were then pooled to purify the PCR products using AMPure XP beads (Beckman Coulter).

Paired-end next generation sequencing

Two biological replicates for Hi-C and capture Hi-C libraries were prepared for each of the cell lines, although sequencing of Hi-C libraries was only conducted on W12 clones with statistically significant virus:host interactions beyond the virus:host breakpoints (i.e. clones D2 and G2). Sequencing was performed on Illumina HiSeq 2500 generating 50 bp paired-end reads (Sequencing Facility, Babraham Institute). CASAVA software (v1.8.2, Illumina) was used to make base calls and reads failing Illumina filters were removed before further analysis. Output FASTQ sequences were mapped to the human reference genome (GRCh37/hg19) containing the HPV16 genome as an extra Chromosome and were filtered to remove experimental artefacts using the Hi-C User Pipeline[82].

HiCUP and SeqMonk

Sequence data was obtained from Illumina HiSeq paired-end sequencing. Using the HiCUP Pipeline[82] paired-end Capture Hi-C (cHi-C) fastq files were mapped with Bowtie 2[83] to a human GRCh37 reference containing a HPV16 pseudo-chromosome. HiCUP removes invalid and artefactual di-tags by overlaying the di-tags on an in silico restriction digest of the reference. The resulting BAM files contained putative di-tags for use in subsequent analyses. SeqMonk (https://www.bioinformatics.babraham.ac.uk/projects/seqmonk/) was used to quantitate and visualise the density of di-tags contained in the BAM files. The HPV16 sequence and annotation files were downloaded from the European Nucleotide Archive (www.ebi.ac.uk/ena/data). ENCODE Annotation for NHEK [84] was obtained from Ensembl release 75 [85].

Circos

The raw cHi-C fastq files were converted to fasta format and BLAST [86] was used to search the HPV16 genome for reads mapping to it. The partner human reads were determined and the 2 sets of reads were mapped to the GRCh37 reference containing the HPV16 pseudo-chromosome using Bowtie 2. The BAM outputs were converted to BED format and modified to be compatible with the circular visualisation tool Circos [87]. The HPV16 genome was split into bins of 500 bp and the count per bin determined from the chimeric human-HPV16 di-tags. The counts, the HPV16 MboI restriction map and gene coordinates were annotated on the Circos plots.

GOTHiC

The HiCUP output was converted to format compatible with the Bioconductor package GOTHiC[88]. To find significant interactions between distal locations GOTHiC implements a cumulative binomial test based on read depth. This was used to identify regions of the human genome in contact with the HPV16 pseudo-chromosome at a resolution of 1kb. Statistically significant interactions were determined by a cumulative binomial test; results were filtered on their adjusted p-value (the q-value), where the threshold in this case was set at <0.05. As such, q-values ranged anywhere between 0.049 to 3.043e-13. Di-tag mappings were visualised with Circos after filtering the previous Circos input by the GOTHiC determined interactions.

Insulation score plots

Measurement of the topological domain structure along the chromosomes was computed with an average insulation score profile at the TAD boundaries. The insulation score is the standardized -log enrichment of contacts between the downstream and upstream 300kb regions (-log (a/(a+b1+b2)) where ‘a’ is the number of contacts between, and ‘b1’ and ‘b2’ the number of contacts within the upstream and downstream 300kb regions. Using this definition, a more positive insulation score indicates a stronger TAD boundary.

RNA-seq analysis and alignment

Complementary DNA (cDNA) libraries were prepared for the five W12 clones (F, A5, D2, H, G2) epigenetically characterised previously[21], two further clones with low HPV16 copy number, distinct integration sites and statistically differing growth rates versus parental lines (R2, 10q22.1, insignificant; B, 8q24.21, significant)[15] to give greater power to statistical analysis here, and the original episomal W12 (Par1) cell line (two biological replicates each) by total RNA extraction from confluent cells with Ribo-Zero rRNA depletion and DNAse treatment before cDNA was prepared with the TruSeq RNA and DNA Sample PrepKit (Illumina). 50bp paired-end cDNA libraries sequenced on an Illumina HiSeq 2000 (Genomics Core Facility, EMBL Heidelberg). Sequence adapters were trimmed from the reads with Kraken[97]. Trimmed FASTQs were mapped against a GRCh37.p13 reference transcriptome (Ensembl version75) that included HPV16 transcript annotation using STAR[98] with default parameters. Strand specific gene counts were obtained from alignments with HTSeq[99] and differential gene expression analysis performed using the R/Bioconductor package DESeq2[100]. Modulation of host transcript levels due to virus integration was then evaluated per clone in comparison to the mean expression of all other clones with RNA-seq libraries (i.e. versus a 6-clone mean expression level), and additionally in comparison to the original episomal W12 cell line. Variance significance (false discovery rate (FDR) of 0.05) was calculated using an ‘F-test of equality of variances’ with a Benjamini Hochberg adjustment of P-values by FDR multiple testing correction.