HPV-CCDC106 integration alters local chromosome architecture and hijacks an enhancer by three-dimensional genome structure remodeling in cervical cancer

1. Introduction

Cervical cancer is the second leading cause of cancer deaths in women in developing countries (Siegel et al., 2017). Human papillomavirus (HPV) infection is the major risk factor, causing 99% of cervical cancers (Walboomers et al., 1999). The process of HPV-associated carcinogenesis is lengthy, with integration of HPV DNA into the host cell genome considered a critical step (Wentzensen et al., 2004). Evidence has shown that HPV integration increases from precancerous cells to cancerous cells. For instance, HPV integration events are found in 53.8% of cervical intraepithelial neoplasias (CINs) and in 81.7% of cervical cancers (Hu et al., 2015b). A classical HPV integration event occurs at human chromosome 8q24, which is located ∼500 kb upstream of the proto-oncogene MYC (Peter et al., 2006; Adey et al., 2013). Using chromosome conformation capture (3C) technology, Shen et al. (2017) showed long-range chromatin interaction between an integrated HPV fragment at 8q24 and the MYC gene. These results highlight the need to study both HPV integration loci and chromatin spatial organization to understand the role of HPV integration in cervical cancer.

HPV integration affects not only the generation of altered transcripts (Cancer Genome Atlas Research et al., 2017) but also interchromosomal rearrangements, leading to genomic instability (Waggoner, 2003; Akagi et al., 2014). Several studies have shown that the molecular pathogenesis of HPV-driven cervical cancer is based on episomal HPV E6/E7 expression and HPV integration (Androphy et al., 1987; Koneva et al., 2018; Oyervides-Munoz et al., 2018). However, previous studies have not validated the role of episomal HPV E6/E7 and HPV integration separately. Furthermore, the comprehensive effects of HPV integration on the host genome structure are not completely understood.

In this study, we collected 61 cervical intraepithelial neoplasia and cancer samples and identified an HPV integration hot spot in the CCDC106 gene using high-throughput viral integration detection (HIVID) sequencing. To investigate the mechanism and consequences of this integration event, we selected fresh cervical cancer tissue containing the unique HPV-CCDC106 integration site at chromosome 19, without nonintegrated episomes, as per the HPV-CCDC106 integration hot spot identified in the study cohort. For the first time, combining linear genome sequencing (i.e., whole-genome sequencing [WGS], RNA sequencing [RNA-seq], and chromatin immunoprecipitation sequencing [ChIP-seq]) and high-throughput chromosome conformation capture (Hi-C), we identified a correlation among chromatin structure, gene expression, and HPV-CCDC106 integrated loci in cervical cancer. These results will help improve our understanding of HPV-CCDC106 integration in cervical carcinogenesis.

2. Results

2.1 HPV integration enrichment on chromosome 19, particularly CCDC106, in cervical cancer

Under all relevant regulations of the Tongji Hospital Institutional Review Board, we collected 61 samples of exfoliated vaginal cells (ThinPrep cytologic test [TCT]) diagnosed with low-grade squamous intraepithelial lesion (LSIL), high-grade squamous intraepithelial lesion, or squamous cell carcinoma (SCC). Using HIVID sequencing (Li et al., 2013; Hu et al., 2015b) to identify HPV integration sites in the 61 samples, we found an HPV-integrated frequency in chromosome 19 of 22.95% (14 samples showed 25 different integration sites, Fig. 1A and Table S1). Of the 14 samples with integration, 10 exhibited HPV integration in the CCDC106 gene (Fig. 1B and Table S1). This HPV integration site is a previously unidentified integration hot spot.

Fig.  1.  HPV integration was enriched on chromosome 19, particularly the CCDC106 gene, in cervical cancer. A: Histogram of HPV16 integration events on chromosome 19 as observed in TCT samples. Cervical precancerous lesions ( n = 40) and cervical cancer (n = 21). B: Distribution of HPV integrations along chromosome 19. Of 14 samples, including a total of 25 integration sites on chromosome 19, 10 displayed HPV integration in CCDC106. C: Circos plot showing HPV integration via HIVID testing for a patient who underwent two cervical biopsies 12 months apart. Right fan diagram is an enlarged representation of chromosome 19. The inner blue circle indicates chronic cervicitis sample result, middle yellow circle indicates HSIL sample results, and outer circle indicates chromosomes. D: FISH results for overlap of CCDC106 (red) and HPV sequences (green) in Carcinoma∗ sample and H&E staining results. HPV-CCDC106 integration is indicated by arrows (of 26 cells counted, 24 showed overlap). E: Quantitative intensity of FISH probe shown with a light blue arrow in D. Red: CCDC106 probe; green: HPV L1 probe; blue: DAPI. X-axis: length of line, Y-axis: fluorescence intensity. HPV, human papillomavirus; TCT, ThinPrep cytologic test; HIVID, high-throughput viral integration detection; HSIL, high-grade squamous intraepithelial lesion; FISH, fluorescence in situ hybridization; H&E, hematoxylin and eosin; DAPI, 4′6-diamidino-2-phenylindole.

DownLoad: Full-Size Img PowerPoint

To verify HPV-CCDC106 integration, we amplified the integration loci from 10 samples using polymerase chain reaction (PCR) and sequenced the DNA using Sanger sequencing (Fig. S1). The sequencing results were consistent with those of HIVID. Interestingly, in two cervical biopsies collected 12 months apart from the same patient with HPV infection, we found the same integration site and same HPV type in both the normal and precancer samples in the CCDC106 gene (Figs. 1C and S2, Table S1), indicating that HPV-CCDC106 integration may persist during cervical cancer progression in high-grade squamous intraepithelial lesion lesions.

To explore the role of HPV-CCDC106 integration in cervical carcinogenesis, we selected a fresh sample containing the HPV-CCDC106 integration site on chromosome 19 (named Carcinoma∗) using HIVID (Table S1) and fluorescence in situ hybridization (FISH) (Fig. 1D and E). The following analyses were based on this clinical sample.

2.2 HPV integration in CCDC106 in carcinoma∗

To explore the effects of HPV-CCDC106 integration on genomic structure, we performed WGS on the Carcinoma∗ sample. By analyzing the WGS data, we found a fragment of HPV16 (coordinates from 3301 to 7906/1–1811, length 6416 bp) integrated on chromosome 19 at genomic coordinates 56161396–56161398 (fourth intron of CCDC106 gene) in the human genome (Fig. 2A and B). Parts of the HPV16 E2 and E1 genes, along with the complete E4, E5, L2, L1, E6, and E7 genes, were integrated into the genome of the Carcinoma∗ sample. A microhomologous “AGT” fragment was present upstream and a microhomologous “CAGC” fragment was present downstream of the HPV integration, similar to the microhomologs reported in our previous study (Hu et al., 2015b).

Fig.  2.  HPV integration in CCDC106 gene in Carcinoma∗ sample. A: Gene model of HPV-CCDC106 integration (yellow and green segments of HPV genome are integrated into Carcinoma∗ sample; gray segment is nonintegrated fragment; orange ‘AGT’ and red ‘CAGC’ indicate microhomologous fragments). B: Chromosome 19 (position 56161396) was linked to HPV 1811 locus, and chromosome 19 (position 56161398) was linked to HPV 3301 locus. C: Sequencing coverage of HPV genome in Carcinoma∗ sample data (red track: HPV sequence coverage; genome data BAM: reads mapped to genome; red indicates reads mapped to forward strand, blue means reverse; blue line indicates split reads or paired-end reads; red gene track: genes of HPV). D: Circos plot of CNVs, indels, and SNPs from Carcinoma∗ sample. Inner circle to outer circle represents CNVs, indels, SNPs, and gene annotations. E: Density of variations on chromosomes; Y-axis represents the number of variations on each chromosome divided by the length of chromosome (orange: indel, purple: SNP). HPV, human papillomavirus; CNVs, copy number variations.

DownLoad: Full-Size Img PowerPoint

To examine whether both HPV integration and extrachromosomal HPV episome replication or only HPV integration occurred in the Carcinoma∗ sample, we examined sequencing depth along the HPV16 genome. No reads were located in the HPV16:1811–3301 region, and the sequencing coverage in the remaining portion of the HPV16 genome containing reads was approximately 8× (Fig. 2C), indicating the absence of extrachromosomal HPV episomes in the Carcinoma∗ sample. We analyzed the reads partially aligned to HPV16:1811 and HPV16:3301 and found that the other sections of these reads were aligned to the human genome, indicating that these reads covered the HPV-human genome integration junctions. These results were confirmed by PCR (Fig. S3 and S4). We used SAMtools to determine sequencing depth near the integration site (Fig. S5). Reads that passed Chr19:56161397 were the nonintegrated chromosome, with a sequencing depth of 26X. Total reads around Chr19:56161397 showed a sequencing depth of 34X, and reads with HPV integration showed a sequencing depth of 8X. Based on read depth, the integration ratio was ∼1:3. Therefore, we concluded that, in the Carcinoma∗ sample, HPV16 integrated into the host genome without extrachromosomal HPV episomes. This high-quality material was thus suitable for studying the effects of pure HPV integration into the human genome, unlike materials with a mix of integrated HPV and extrachromosomal HPV episomes.

We then focused on whether HPV integration could induce instability in the entire genomic structure, which can produce substantial variation (Zhao et al., 2016). We identified 99136 single-nucleotide polymorphisms (SNPs), 369168 insertions or deletions (indels), and 8966 copy number variations in the WGS data (Fig. 2D and Tables S2–4). Via annotations in the Single-Nucleotide Polymorphism Database (dbSNP) (Sherry et al., 2001) and the Catalogue of Somatic Mutations in Cancer (COSMIC) database (Forbes et al., 2010, 2011), we identified 12 SNPs associated with cervical cancer (i.e., COSMIC labels: COSM1134142, COSM111734, COSM1044560, COSM1044558, COSM1595268, COSM1088016, COSM1083212, COSM1598410, COSM943980, COSM1477981, COSM1477980, and COSM1588240). Interestingly, we observed a higher density of mutation genes (SNPs and indels per Mb) on chromosome 19 (Fig. 2E), which was normalized by the length of each chromosome.

2.3 Transcriptional changes in HPV-CCDC106 integration in Carcinoma∗ cells compared with primary cultured normal cervical epithelial cells

We next examined the transcription levels of RNA in Carcinoma∗ and normal cervical epithelial (Ncepi) cells. Ncepi cells are primary cultured cells derived from the normal cervix (see Materials and Methods). Based on RNA-Seq data, we identified differentially expressed genes (DEGs) between Carcinoma∗ and Ncepi cells using differential expression analysis (Fig. 3A and Table S5). We found that the expression levels of genes in the tumor suppressor BRCA family and CAV family, as well as PEG3 (near the integration site), were significantly decreased. Gene Ontology (GO) (Fig. S6) and KEGG pathway analyses (Fig. S7 and Tables S6–7) of the DEGs identified several significantly enriched pathways, including positive regulation of leukocyte cell-cell adhesion and the p53 signaling pathway.

Fig.  3.  HPV-CCDC106 integration affected transcription of many cancer-related genes. A: Differentially expressed genes (DEGs); blue dots indicate downregulated genes, and red dots indicate upregulated genes (CCDC106# is the sum of CCDC106 gene and CCDC106 fusion gene). B: Schematic of fusion gene formed by HPV and CCDC106. Sanger sequencing around junction points is shown at bottom. C: RNA sequencing coverage for Carcinoma∗ and Ncepi cells. Y-axis represents RNA coverage (coverageˆ0.3). D: Comparison of expression level of CCDC106 estimated as covering fragments from exons 4 and 5 in Carcinoma∗ and Ncepi cells (counts of split reads or paired-end reads from exon 4 to exon 5 of nonfusedCCDC106 gene, normalized by the total number of reads). E: Real-time PCR analysis of expression levels of exons 4–5, Fusion-1, and Fusion-2 in Ncepi and Carcinoma∗ cells (N/A, not available. Each group included two replicates, represented as mean ± SD). F: Percentage of DEGs on each chromosome; Y-axis represents percentage of upregulated (red) and downregulated genes (blue) on each chromosome. G and H: Distribution of DEGs on chromosome 19 (G) and near integration site (H) (red: upregulated genes, blue: downregulated genes, green: HPV integration site). HPV, human papillomavirus; PCR, polymerase chain reaction; Ncepi, normal cervical epithelial; SD, standard deviation.

DownLoad: Full-Size Img PowerPoint

RNA-Seq indicated that the expression of CCDC106 (i.e., total expression of different types of CCDC106, including CCDC106 fusion genes) increased significantly (Fig. 3A). Based on STAR-Fusion (version 1.5.0) (Akers et al., 2018), we identified 31 fusion genes from the RNA-Seq data (Table S8). Of note, we identified fusion transcripts of CCDC106 and HPV16 (Fig. 3B and Methods) and spliced these fusion genes using StringTie (version 1.3.3) (Pertea et al., 2016) and Trinity (version 2.8.4) (Grabherr et al., 2011) (see Materials and Methods). One fusion gene included the first four exons of CCDC106 and a portion of the HPV16 E2 gene; the other included the E6, E7, and E1 genes of HPV16 and the fifth exon of CCDC106, which showed three alternative splicing transcripts (Fig. S8). These fusion genes were validated by PCR and Sanger sequencing (Fig. 3B). In general, we identified the nonfused CCDC106 gene (junction of exon 4 and 5 of CCDC106) and HPV-CCDC106 fusion gene in the RNA-Seq data (Figs. 3B, 3C and S8). We then compared the expression levels of the nonfused CCDC106 gene in the Carcinoma∗ and Ncepi samples (Fig. 3D). Results demonstrated lower expression levels of the nonfused CCDC106 gene in the Carcinoma∗ sample, as confirmed by real-time PCR (Fig. 3E).

Moreover, the density of DEGs, especially upregulated DEGs, was higher on chromosome 19 than on any other chromosome (Fig. 3F). We also found that the number of DEGs increased closer to the integration site (Fig. 3G–H and Tables S9–10).

2.4 HPV-CCDC106 integration altered local 3D genome structure in Carcinoma∗ cells

To demonstrate the effects of HPV-CCDC106 integration on the 3D structure of the host genome, we applied our newly developed digestion-ligation-only Hi-C (DLO Hi-C) method (Lin et al., 2018) to the Ncepi and Carcinoma∗ cells. The overall heatmap revealed higher order genome organization (Fig. 4A). The Carcinoma∗ cell heatmap was similar to that of the Ncepi cells, except for several obvious abnormal interactions, which were possibly caused by chromatin translocation (Fig. 4A). The long-range interaction signal of the Carcinoma∗ cells was stronger than that of the Ncepi cells (Table S11). Interactions among checkerboard signals inside chromatin on chromosome 19 can represent open and closed compartments (Lieberman-Aiden et al., 2009). Here, compared with the signals from the Ncepi cells, the checkerboard signals from the Carcinoma∗ cells were enhanced and several differences in the compartmentalization pattern could be observed (Fig. 4B). We found that most A/B compartments remained conserved between the Ncepi and Carcinoma∗ cells, with only 21% of the annotated genome changing from A to B or B to A compartments (Fig. S9). Together with the gene expression data, these results indicate that the expression levels of genes changing from B to A tended to increase, whereas those changing from A to B tended to decrease (Fig. S10 and Tables S12–14), consistent with previous research (Sade-Feldman et al., 2017). Interestingly, based on the number of DEGs divided by the number of genes changing compartments on each chromosome, chromosome 19 had the highest proportion of genes changing from A to B and chromosomes 19 and 21 had the highest proportion of genes changing from B to A (Fig. 4C). However, only a small number of genes (nine) changed from B to A on chromosome 21 (Tables S5 and S12), potentially contributing to the higher proportion value for this chromosome.

Fig.  4.  HPV-CCDC106 integration altered local 3D genome structure of Ncepi and Carcinoma∗ cells. A: Heatmaps of whole genome (resolution: 1 Mb) and chromosome 19 (resolution: 100 kb) in Ncepi and Carcinoma∗ cells. B: Heatmaps indicate contact correlation matrix for chromosome 19 in Ncepi and Carcinoma∗ cells. Corresponding A/B compartments are shown at bottom. C: Bubble diagram showing percentage of differentially expressed genes (DEGs) moving between compartments. Sizes of bubbles indicate the number of genes. D: Pie chart showing percentages of unchanged TADs and changed TADs in Ncepi and Carcinoma∗ cells. Bar plots of the number of changed TADs (E) and the number of DEGs (F) in changed TADs per Mb. Circos plot of chromosome 19 (G) and enlarged region (Chr19:40–58 M) (H) (track from inside to outside represents expression levels in Carcinoma∗ and Ncepi∗ cells, distribution of SNPs, indels, CNVs in Carcinoma∗ cells, gene density, and cytobands. Links: differential interactions of integration region [Chr19:55–57 M; red triangle]).I: Circos plot of immune-related genes (outside text) that interacted differentially with integration region (links: intrachromosomal interactions at minimum distance of 500 kb between corresponding regions [resolution, 100 kb] with a frequency of >8).J: Heatmap of Z-score-transformed FPKM values. Depicted genes are those differentially interacting with integration region (Chr19:55–57 M). HPV, human papillomavirus; FPKM, fragments per kilobase of exon model per million mapped reads; TADs, topologically associated domains; CNVs, copy number variations; SNPs, single-nucleotide polymorphisms.

DownLoad: Full-Size Img PowerPoint

Topologically associated domains (TADs) are important stable regulatory units closely related to biological processes such as replication (Pope et al., 2014). The number and average size of TADs did not differ significantly between Ncepi and Carcinoma∗ cells (Fig. S11). However, several local areas exhibited significant differences. Among the identified TADs (approximately 4000), ∼55% were unchanged, but ∼45% demonstrated some form of modification, including fusion, separation, and shift (Fig. 4D). Among all chromosomes, chromosome 19 had the highest number of changed TADs per Mb (Fig. 4E; Table S15). In addition, both the number of DEGs per Mb and the percentage of DEGs within the boundaries of changed TADs across all genes were appreciably higher on chromosome 19 than on any other chromosome (Figs. 4F and S12).

By comparing the 3D structures of Ncepi and Carcinoma∗ cells, we found that some genomic regions exhibited differential interactions on chromosome 19 (Fig. 4G and Table S16), especially near the integration site (Fig. 4H). By comparing the genomic regions that interacted with the HPV-CCDC106 integration site, we found many immune-related genes that showed differential expression between the two samples (Fig. 4I and J).

2.5 Enhancer hijacking by HPV-CCDC106 integration downregulated tumor suppressor gene PEG3 but upregulated CCDC106 via TAD structure remodeling

We next focused on the 3D structure around the HPV-CCDC106 integration site. The integration site inCCDC106 was located at the boundary between two TADs. The left TAD was the same in both Ncepi and Carcinoma∗ cells, whereas the right TAD in the Ncepi cells was divided into two smaller TADs in the Carcinoma∗ cells (Fig. 5A). A conserved enhancer (Fig. 5C) annotated in NHEK and HeLa–S3 cells significantly interacted with the PEG3 promoter in the Ncepi cells (Fig. 5D). After TAD division in the Carcinoma∗ cells, the interaction between the enhancer and PEG3 gene promoter was significantly reduced, whereas the interaction between the enhancer and CCDC106 gene promoter was significantly increased (Fig. 5D). Therefore, we hypothesized the following model (Fig. 5B): 1) HPV-CCDC106 integration hijacked an enhancer via division of the TAD structure; 2) after TAD division, the new TAD boundary acted as an insulator to block enhancer regulation of PEG3, resulting in its downregulation; 3) the enhancer then interacted with the CCDC106 promoter and upregulated the expression of CCDC106. This model was validated by virtual 4C analysis of Hi-C data and epigenetic patterns (Fig. 5D). Subsequently, the protein expression levels of PEG3 and CCDC106 were assessed by immunohistochemical (IHC) staining in the Carcinoma∗ tissue and in adjacent cervical tissue. Compared with adjacent cervical cells, the Carcinoma∗ cells exhibited upregulated CCDC106 expression and downregulated PEG3 expression. These expression patterns were validated using 10 biopsy samples containing the HPV-CCDC106 integration (validated in Figs. 1B and S1, Table S1) and 10 biopsy samples without the HPV-CCDC106 integration from our cohort (Figs. 5F and S13). Results demonstrated that the expression levels of CCDC106 and PEG3 were upregulated and downregulated, respectively, in samples showing the same HPV-CCDC106 integration.

Fig.  5.  Enhancer hijacking around integration site by HPV-CCDC106 integration via TAD structure remodeling. A: Hi-C chromatin interactions in Chr19:55.52M–57.68M region of Ncepi and Carcinoma∗ cells (resolution: 40 kb) (black triangle: TAD regions; red arrow: HPV-CCDC106 integration site). B: Model illustrating change in chromatin TADs caused by HPV-CCDC106 integration (triangle: TAD region; blue arrow: PEG3 gene; yellow ellipse: enhancer; blue hexagon: TAD boundary; green arc arrow: chromatin interaction). C: Illustration of epigenetic marks in Chr19:56.39M–58.88M region in NHEK and HeLa–S3 cell lines. ChIP-Seq data were obtained from ENCODE. D: Heatmaps of differential interactions of Chr19:55M–58M region in Ncepi and Carcinoma∗ cells. Middle rectangle indicates putative enhancer. Bottom histogram shows pseudo-4C profiles from DLO Hi-C analysis for putative enhancer region (gray: coverage of normalized fragments; red: regions that significantly interact with viewpoints; blue: background distribution). E: IHC staining of Carcinoma∗ and adjacent cervical samples. F: IHC staining analysis of 10 biopsy samples with HPV-CCDC106 integration and 10 biopsy samples without HPV-CCDC106 integration. ∗∗P < 0.01, ∗∗∗P < 0.001. Student'st-test. Means and SD are shown. HPV, human papillomavirus; SD, standard deviation; IHC, immunohistochemical; TAD, topologically associated domain; Ncepi, normal cervical epithelial; DLO, digestion-ligation-only.

DownLoad: Full-Size Img PowerPoint

To study the interaction patterns around the TAD boundaries, we checked the CTCF binding sites obtained from ENCODE data. After checking the CTCF binding sites of the NHEK cell line, we found the CTCF was bound to the nearest upstream boundary region (near chr19:56900000) (Fig. S14). In the Ncepi cells, there may be CTCF binding sites in the region (near chr19:56900000). However, based on Hi-C data analysis, the region near chr19:56900000 was not a boundary region in the Ncepi cells. By analyzing the HPV-CCDC106 integration region, we found that CCDC106 could combine with the specific enhancer due to development of TAD boundary barriers after new TAD formation.

After TAD formation, we also found that the expression levels of genes inside TAD1 in the Carcinoma∗ cells were generally higher than those in the normal cervical cells. The expression levels of most genes inside TAD2 decreased. The expression levels of genes inside TAD3 tended to remain unchanged or decreased (Fig. S15). In the normal cells, an enhancer acted on the PEG3 gene but acted on CCDC106 after HPV-CCDC106 integration via remodeling of the 3D genome conformation. In addition, after HPV integration, the hijacking of the enhancer from TAD2 and TAD3 to TAD1 contributed to the expression dysregulation observed in the Carcinoma cells.

3. Discussion

HPV integration is the main cause of cervical cancer. Based on all previously reported HPV integration sites (Hu et al., 2015b; Bodelon et al., 2016), we calculated the average HPV integration density per 100 Mb on each chromosome and found that chromosome 19 had the highest incidence of HPV integration events (Fig. S16). The CCDC106 gene on chromosome 19 was then identified as an integration hot spot in our TCT cohort. These findings suggest that HPV-CCDC106 integration may play a crucial role in cervical carcinogenesis, although additional research is required.

For further analysis, we used a fresh cervical cancer tissue sample that showed the unique HPV-CCDC106 integration but no extrachromosomal HPV episomes. Previous research has reported that HPV integrations result in higher expression of genes surrounding the integration loci, as evidenced by host-virus fusion (Khoury et al., 2013; Cancer Genome Atlas Research et al., 2017) and transcriptome fusion (Wu et al., 2018). Here, we not only observed an increase in the integrated gene but also identified enrichment of genome variation and expression change densities in the HPV-CCDC106 integrated chromosome. The genes located near the integration site exhibited more obvious changes than those further away. Furthermore, we identified many upregulated and downregulated genes compared with normal cervical cells. Indeed, CCDC106 and PEG3 expression levels were significantly changed in Carcinoma∗ cells compared with normal cervical cells, consistent with previous findings that HPV integration elevates gene expression at the integration site (Tang et al., 2013; Ojesina et al., 2014). However, these changes cannot be explained simply by the linear genome sequence (Politz et al., 2003).

As the human genome exhibits 3D organization within the cell (Cremer and Cremer, 2001), the complex relationships among chromatin structure, HPV integration functional state, and gene activity require 3D sequencing techniques, including 3C or one of its variants, such as Hi-C (Lieberman-Aiden et al., 2009), or chromatin interaction analysis with paired-end tag sequencing (ChIA-PET) (Fullwood and Ruan, 2009). Here, we obtained Hi-C data from the Carcinoma∗ and Ncepi cells, which showed that gene expression was affected by the 3D structure, including A/B compartments, TADs, and interactions. In particular, compared with that on other chromosomes, TAD modifications per Mb were more obvious on chromosome 19, which also featured HPV integration, consistent with the enrichment of DEGs on chromosome 19. Furthermore, by studying the effects of 3D remodeling on the genomic region near the HPV16-CCDC106 integration site, we observed that the stability of the local TAD was affected, causing TAD shuffling (Spielmann et al., 2018) after HPV sequence integration into the TAD boundary. The downstream TAD was also separated into two smaller TADs. These structural variations across TAD boundaries may have contributed to the fusion of two regulatory domains that do not naturally belong together, causing enhancer adoption. Furthermore, these variations could lead to regulatory loss of function by removing the enhancer from its target genes. Several previous studies using CRISPR/Cas9 to edit DNA sequences have shown changes in 3D genome structures, including enhancer-promoter interactions, which can reach 100 kb in genomic length (Guo et al., 2015; Lupianez et al., 2015).

PEG3 shows tumor suppressor activity in many types of cancer ( Kohda et al., 2001; Feng et al., 2008) and contributes to p53 activation (Relaix et al., 2000). CCDC106 is a p53-interacting partner in vivo (Vogelstein et al., 2000) and is reported to promote the degradation of the p53 protein in cervical cancer cell lines (Zhou et al., 2010). Furthermore, CCDC106 knockdown or inhibition of CCDC106 phosphorylation can suppress tumor progression in breast and cervical cancer (Ning et al., 2019). Here, the Carcinoma∗ sample only contained unique integration loci at CCDC106 (19q13.42, Chr19/hg19:56161397) with no HPV episomal DNA. The 19q13.42 integration of HPV in cervical cancer has been detected in previous studies. Specifically, Yu et al. (2005) and Hu et al. (2015b) reported one (out of 36) and four (out of 103) cases of 19q13.42 integration, respectively, with the same integration site in CCDC106 (Chr19/hg19:56161488) (Fig. S17A). We also found high expression of CCDC106 to be associated with poor disease-free survival (DFS) in cervical cancer (p (log-rank) = 0.0077, hazard ratio (HR high) = 2.2, p (HR) = 0.0094) based on The Cancer Genome Atlas (TCGA) data set (Fig. S17B). These findings indicate that HPV-CCDC106 integration may influence genomic accessibility on integrated chromosomes and dysregulate the expression of nearby genes in carcinogenesis.

Despite the positive results, our research has several limitations. HPV integration can dysregulate the expression of oncogenes and promote carcinogenesis in several ways (McBride and Warburton, 2017), among which HPV enhancer hijacking is only one. Our new method, which considered 3D architecture, only showed that HPV DNA hijacked an adjacent enhancer to drive dysregulation of host genes. In addition, the number of HPV integration sites can vary among samples (Hu et al., 2015b), even if CCDC106 is a hot spot. Thus, our findings based on HPV-CCDC106 integration only partially account for the ways in which HPV may promote carcinogenesis. In addition, the common mechanism that HPV integration alters gene expression and hijacks an enhancer by 3D structural remodeling needs subsequent Hi-C analyses covering additional HPV integration sites.

In summary, we identified a new, previously unreported HPV integration site (CCDC106) and demonstrated enrichment of genome variation and gene expression change in the HPV-CCDC106 integrated chromosome of a cervical carcinoma sample. We observed a correlation between gene expression and 3D structural change of chromosomes, suggesting that HPV-CCDC106 integration altered the 3D genome structure and gene expression in the associated chromosome. In addition, HPV-CCDC106 integration led to the separation of TAD structures, which changed the interactions between the enhancer and host genes. These findings shed new light on the importance of the 3D human genome structure when studying the effects of HPV integration in cervical carcinogenesis.

4. Materials and methods

4.1 Sample preparation and description

All tumor specimens, TCT samples, and normal cervical tissues were obtained from the Department of Gynecological Oncology at Tongji Hospital (Tongji Medical College, Huazhong University of Science and Technology, Wuhan, Hubei, China) under ethics board approval (TJ-IRB20180611). Documented informed consent was obtained from all patients. Surgically resected specimens or TCT samples were immediately frozen in liquid nitrogen and stored at −80 °C for less than one year. DNA and mRNA were extracted from samples that contained a >50% malignant epithelial cell component, as validated by hematoxylin and eosin staining. A detailed description of the TCT samples is available in Table S17.

4.2 Primary cervical epithelial cell (Ncepi) culture and identification

Primary derivative culture of Ncepi cells was performed as described previously (Hines and Allen-Hoffmann, 1996; Bononi et al., 2012). Briefly, we collected surgically resected cervices from 10 women diagnosed with uterine prolapse and treated with hysterectomy. The epithelium was first dissociated using 2.5 U/mL dispase II (04942078001; Roche, Switzerland) at 4 °C overnight and then dispersed in 0.25% trypsin (12605010 Gibco, Denmark). The epithelial cells were counted and cultured in keratinocyte serum-free medium (KSFM, 10785-012; Invitrogen, USA) supplemented with 2.5 μg of human recombinant epidermal growth factor (PHG0311; Gibco, USA) and antibiotics (penicillin 100 U/mL, streptomycin 100 μg/mL, 10015140163; Gibco) at 37 °C with 5% CO₂. All cells used for sequencing were tested for HPV negativity using a Cobas 4800 system (Roche) at the Department of Gynecological Oncology at Tongji Hospital and were identified by fluorescence-activated cell sorting (FACS) to verify an epithelial cell component of >94% (CD326, 324204; BioLegend, USA). All identification testing results are available in Fig. S18.

4.3 Annotation of HPV integration sites

To identify breakpoints in the human genome and HPV integration sites, we carried out HIVID in accordance with previously described algorithms and protocols (Li et al., 2013; Hu et al., 2015b). Detailed data regarding the HPV integration sites in all samples are available in Table S3.

4.4 Validation of HPV integration sites

PCR and Sanger sequencing were carried out to verify the HPV integration sites identified by HIVID. The PCR primers were designed with Primer3 (http://bioinfo.ut.ee/primer3-0.4.0/primer3/); one primer was in CCDC106, and the other was located in the HPV16 genome (Fig. S4). All PCR protocols were performed in accordance with the manufacturer's instructions using a GeneAmp PCR System 9700 thermal cycler (Bio-Rad, USA). The resulting DNA was sequenced on an Applied Biosystems DNA analyzer (Life, Inc., USA).

4.5 Identification of control cell line

Human embryonic lung fibroblasts (MRC5), an identified control cell line for Ncepi cells, were obtained from the Cell Bank of the Chinese Academy of Sciences. The MRC5 cells were cultured in minimum Eagle's medium (MEM) medium supplemented with 10% fetal bovine serum (10099-141; Gibco) and antibiotics (penicillin-streptomycin, 15140163; Gibco). Cells were grown at 37 °C with 5% CO₂.

4.6 Immunofluorescence

Immunofluorescence was performed as described in earlier research (Ding et al., 2014). Ncepi and MRC5 cells (control cell line) were grown on glass coverslips and then fixed, permeabilized, incubated with antibodies (37 °C for 1 h), and imaged by fluorescence microscopy for identification. The antibodies included anti-pan-keratin (4545T; CST, USA), anti-vimentin (5741T; CST, USA), and anti-fluorophore (488/594; Antgene, China).

4.7 Fluorescence in situ hybridization

FISH was carried out as described previously (Ding et al., 2014; Hu et al., 2015a). Briefly, paraffin-embedded sample slides were deparaffinized, pretreated with sodium thiocyanate (NaSCN), digested with pepsin, and hybridized with a Cy3-labeled HPV16-L1 probe and fluorescein isothiocyanate (FITC)-labeled CCDC106 probe. Nuclei were stained with 4′6-diamidino-2-phenylindole. Digital images were obtained with a fluorescence microscope (Olympus BX53, Japan) equipped with Cy3, FITC, and 4′6-diamidino-2-phenylindole filters.

4.8 IHC staining

IHC staining was performed as reported previously (Hu et al., 2015a). Formalin-fixed and paraffin-embedded cervical cancer and normal cervical tissues were subjected to staining using an appropriate kit (ZSGB-Bio, Beijing, China) in accordance with the manufacturer's instructions. Antibodies against CCDC106 (abs125584; Absin, China) and PEG3 (GTX37357; GeneTex, USA) were used. The expression levels of PEG3 and CCDC106 protein in the samples were scored by intensity and extent. Staining intensity was scored as negative, 0; weak, 1; moderate, 2; and strong, 3. Staining extent was scored as per the percentage of stained cells in each field (100 cells/field): none, 0; <25%, 1; 25%–50%, 2; 50%–75%, 3; and >75%, 4. The overall staining score for each sample was determined by multiplying the staining intensity by staining extent. To evaluate the expression of PEG3 and CCDC106, lesions were evaluated by high-power (×400) microscopy, with each point in the figure representing the mean score from two pathologists.

4.9 RNA-seq

Total RNA was extracted using RNAiso Plus (9109; Takara, Japan) in accordance with the manufacturer's instructions. Sequencing libraries were prepared using a VAHTS Stranded mRNA-Seq Library Prep Kit (NR602-02; Vazyme, China) as per the manufacturer's protocols.

4.10 Digestion-ligation-only Hi-C

DLO Hi-C was performed as described in earlier research (Lin et al., 2018). First, 5 × 10⁶ cells were cross-linked with 1.5 mmol/L ethylene glycol bis (succinimidyl succinate, EGS, 21565, Thermo, USA) and 1% formaldehyde solution (Sigma, Darmstadt, Germany) and lysed with lysis buffer (10 mmol/L Tris-HCl [pH: 8.0], 10 mmol/L NaCl, 0.3% Igepal CA-630, and complete protease inhibitor [Roche]). The nuclei were then gently resuspended in 100 μL of 1% sodium dodecyl sulphate (SDS) solution and incubated at 60 °C for 5 min. The nuclei were divided into two 1.5-mL tubes (tubes A and B) before the addition of 15 μL ofHindIII (100 U/μL, NEB), 50 μL of Triton X-100 (20%), and 385 μL of NE buffer 2.1, followed by digestion of chromatin at 37 °C for 6 h. After digestion with restriction enzymes, two different half-linkers (i.e., linker A and linker B) were ligated into the digested chromatin. Next, 50 μL of NE buffer 2.1 (10×), 200 μL of half-linkers (linker A and B) (600 ng/μL), 10 μL of 100 mmol/L ATP (Thermo), 10 μL of HindIII (100 U/μL, NEB), 5 μL of T7 DNA ligase (3000 units/μL, NEB), and 225 μL of dd-H₂O were added to 500 μL of digested chromatin in each tube (A and B), followed by incubation at 25 °C and 15 rpm for 1 h. After ligation, the chromatin-protein complex was centrifuged at 17000 rpm for 1 h at 4 °C to remove excess unbound linkers. The supernatant was removed, and the pellet was resuspended in 1 mL of phosphate-buffered saline (PBS) and centrifuged at 17000 rpm for 1 h at 4 °C. The pellets from tubes A and B were then dissolved in 50 μL of 0.5% SDS and combined into one 15-mL tube. Then, 10 mL of preheated (37 °C) T4 DNA ligation buffer (1×, Thermo) containing 0.4% low-melting point agarose gel (Takara), 1% Triton X-100, and 400 units of T4 ligase (EL0011; Thermo) was added to the tube. The gel was mixed, placed in ice water to solidify, and treated at 20 °C for 6 h for proximity ligation. The agarose gel was melted and digested with 20 μL of agarose (Takara) at 65 °C for 1 h and then concentrated to 800 μL using an Amicon-Ultra-0.5 centrifugal filter unit (NMWL, 30 kDa, Millipore, USA). After this, 100 μL of 10% SDS and 100 μL of 10 mg/mL proteinase K (Sigma) were added, followed by incubation for 3 h at 60 °C to fully release DNA. After extracting DNA using standard DNA extraction protocols and dissolving in 160 μL of dd-H₂O, 20 μL of CutSmart buffer (10×), 10 μL of S-adenosylmethionine (SAM, NEB), and 10 μL of MmeI (2 units/μL, NEB) were added to 160 μL of the DNA sample, with incubation performed at 37 °C for 1 h for digestion. An 80-bp contact DNA fragment was then recycled on the agarose gel and dissolved in 40 μL of dd-H₂O. Finally, 3 μL of in house–generated Illumina sequence adapter, 5 μL of T4 DNA ligase buffer (10×), and 2 μL of T4 ligase (2 units/μL, Thermo) were added to 40 μL of the DLO Hi-C DNA fragment and incubated at 16 °C for 30 min. The ligation mixtures were then purified using AM Pure XP beads and eluted in 40 μL of dd-H₂O. The eluted DNA was repaired in PreCR Repair Mix (NEB) at 37 °C for 20 min in a final volume of 50 μL. Using 5–10 μL of the repaired DNA as a template and amplifying for less than 13 cycles, the PCR amplification product (200 bp) was used for sequencing.

4.11 DLO Hi-C pipeline

The DLO Hi-C sequence reads were processed via a previously described pipeline ( Lin et al., 2018). A modified Java program, DLO Hi-C Tool (Hong et al., 2020), was used to identify linker composition (AA, BB, AB, or BA). Only sequences with homodimeric AA or BB linkers were trimmed for further analysis. The paired-end tag (PET) sequences were separately mapped to the human reference genome (hg19) by Burrows-Wheeler Alignment (BWA) (Li and Durbin, 2009). Uniquely mapped reads (mapping quality score [MAPQ]≥20) were paired by BEDtools (Quinlan, 2014). The redundant reads arising from clonal PCR amplification were removed. The nonredundant mapped reads were then classified as either valid or invalid interactions. Self-ligation (reads mapped to the same fragment) and religation (reads mapped to an adjacent fragment with a nearby position) may lead to invalid interactions. The iterative correction (ICE) method (Imakaev et al., 2012) was applied to normalize the interaction matrix at different resolutions. For compartment analysis, Juicer (Durand et al., 2016) was applied to calculate the Pearson correlation matrix at a resolution of 400 kb. For heatmap analysis, before converting the interaction matrix into differential heatmaps, the DLO Hi-C matrix was normalized by sequencing depth. The R package diffHiC (Lun and Smyth, 2015) was used to calculate differential interactions at a resolution of 500 kb.

4.12 Significant interaction calling

To estimate confident intrachromosomal contacts, the matrix was converted to another interaction matrix with a window size of 10 HindIII restriction fragments. Statistical confidence was estimated by Fit-Hi-C (Ay et al., 2014). Interactions with a q-value of ≤0.001 were considered significant.

4.13 Identification of TADs and differential TADs

To identify TADs, interaction matrices were binned at a 40-kb resolution and iteratively corrected using a previously described pipeline (Imakaev et al., 2012). Rscript TopDom (Shin et al., 2016) was used with default parameters to calculate TAD regions. If the boundaries of two TADs were within 80 kb of each other, they were considered unchanged; otherwise, they were defined as shifted. TADs that separated into two smaller TADs were defined as separated. TADs that merged into a larger TAD were defined as fused.

4.14 Calculation of TAD changes on each chromosome

The average number of changed TADs per Mb was calculated as the number of changed TADs divided by the length of each chromosome. The average number of DEGs in changed TADs per Mb was calculated by the number of DEGs in changed TADs divided by the length of each chromosome. The percentage of DEGs in changed TADs was calculated by the number of genes of changed TADs divided by the number of genes of each chromosome.

4.15 Virtual 4C signal generation

Virtual 4C data were generated from the Hi-C data with the viewpoint of Chr19:56580354–56608205. Peaks were called in the virtual 4C data set using peakC, as described previously (Geeven et al., 2018), with the following parameters: minDist = 40000,wSize = 11, qWd = 1.5, and qWr = 1.5.

4.16 WGS analysis

Quality control was performed using FastQC (version 0.11.5) (Brown et al., 2017) and Trimmomatic (version 0.36) (Bolger et al., 2014). We created a new genome file that contained the human genome (hg19 with chromosome Y removed) and HPV16 genome (https://pave.niaid.nih.gov/#explore/reference_genomes/human_genomes/locus_view/fetch?id=HPV16REF&format=Locus%20view&hasStructure=none), with the sequence data then aligned to a new reference genome by SpeedSeq (version 0.1.0) (Chiang et al., 2015) under default parameters. Structural variations were then identified using the SpeedSeq sv module. Local realignment was carried out with the Genome Analysis Toolkit (GATK) (3.8-0-ge9d806836), and SNPs and indels were called with the GATK HaplotypeCaller. The identified SNPs and indels were annotated with ANNOVAR (http://annovar.openbioinformatics.org/en/latest/) (Wang et al., 2010). The copy number variations were called by the SpeedSeq sv module (Chiang et al., 2015). Statistical analysis of data is available in Table S18.

4.17 RNA-seq analysis

The RNA-Seq data underwent quality control with FastQC (version 0.11.5) (Brown et al., 2017) and Trimmomatic (version 0.36) (Bolger et al., 2014). We created a new genome file that contained the human genome (hg19 with chromosome Y removed) and HPV16 genome and aligned the data to the new reference genome using hisat2 (version 2.0.5) (Kim et al., 2015; Pertea et al., 2016) with default parameters. Gene expression was then computed for each gene using HTSeq (version 0.9.1) (Anders et al., 2015), and counts were then transformed to fragments per kilobase of exon model per million mapped reads. DEGs were analyzed by DESeq2 (version 1.14.1) (Love et al., 2014). STAR-Fusion (version 1.5.0) (Akers et al., 2018) was used to identify fusion transcripts. SAMtools (version 1.4.1) (Li and Durbin, 2009) was used to obtain reads mapped to CCDC106 and HPV, with Trinity (version 2.8.4) (Grabherr et al., 2011) used to assemble the reads. We combined the CCDC106 and HPV genome sequences and annotations and used StringTie (version 1.3.3) (Pertea et al., 2016) to find alternative splicing transcripts. Statistical analysis of paired-end RNA-Seq data is available in Table S4, and Spearman correlation of read counts is available in Fig. S19. We used a normal cervical tissue mixed culture as the control. The identification process was developed RNA-Seq data map to hg19 genome with hisat2 (Kim et al., 2015; Pertea et al., 2016), calculated the gene read count by HTSeq (version 0.9.1) (Anders et al., 2015), and identified DEGs by DESeq2 (version 1.14.1) (Love et al., 2014).

4.18 Data visualization

The genomes of samples and HPV from the sequencing data were constructed using Integrative Genomics Viewer (version 2.5.0) (Robinson et al., 2011). The ChIP-Seq result from ENCODE was viewed using the WashU Epigenome Browser (Zhou et al., 2011).

4.19 Primer sequences for PCR and probe sequences for FISH

All sequences of primers or probes are available in Table S20.

4.20 Statistical analysis

Statistical analysis was performed by SPSS, version 12.0. Data are expressed as mean ± SD. A P-value of less than 0.05 indicated a significant difference between the two groups.

4.21 Ethics approval and consent to participate

All tumor specimens, TCT samples, and normal cervical tissues were obtained from the Department of Gynecological Oncology at Tongji Hospital under ethics board approval (TJ-IRB20180611) and with documented informed consent from all patients.

4.22 Data deposition and material sharing

The raw sequencing data supporting the findings of this research were deposited in the Genome Sequence Archive (GSA) under accession number CRA001401 (http://bigd.big.ac.cn/gsa/s/P37lFNi0).

CRediT authorship contribution statement

Canhui Cao: Writing- Original draft preparation, Methodology, Conceptualization, Visualization, Validation, Investigation. Ping Hong: Writing- Original draft preparation, Methodology, Software, Formal analysis. Xingyu Huang: Writing- Original draft preparation, Methodology, Software, Formal analysis. Da Lin: Methodology. Gang Cao: Validation. Liming Wang: Data curation. Bei Feng: Data curation. Ping Wu: Data curation. Hui Shen: Data curation. Qian Xu: Software, Validation. Ci Ren: Data curation. Yifan Meng: Data curation. Wenhua Zhi: Data curation. Ruidi Yu: Data curation. Juncheng Wei: Funding acquisition, Resources. Wencheng Ding: Resources. Xun Tian: Resources. Qinghua Zhang: Resources. Wei Li: Resources., Qinglei Gao: Supervision, Resources, Funding acquisition. Xiaoyuan Huang: Resources. Gang Chen: Resources. Ling Xi: Resources. Wing-Kin Sung: Resources, Writing- Reviewing and Editing. Zheng Hu: Project administration, Supervision, Resources. Hui Wang: Project administration, Supervision, Funding acquisition. Guoliang Li: Supervision, Project administration, Funding acquisition, Writing- Reviewing and Editing. Peng Wu: Project administration, Funding acquisition, Writing- Reviewing and Editing.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (81630060 to P.W., 31771402 to G.L., 81830074 and 81772786 to H.W., 81572569 to G.C., and 81772775 to J.W.), National Science and Technology Major Project (2019YFC1005202 and 2019YFC1005201 to K.L., and 2018ZX10301402-002 to Q.G.), and the research-oriented clinician funding program of Tongji Medical College, Huazhong University of Science and Technology for P.W. We also thank the ENCODE for public available data.