Human cervical cell lines SiHa (ATCC HTB-35) and Ect1/E6E7 (ATCC CRL-2614) were obtained from the American Type Culture Collection (ATCC). Normal human epidermal keratinocytes (NHEK) were purchased from PromoCell (C-123003; Germany). All the cells were authenticated by short tandem repeat (STR) tests. SiHa/NHEK and Ect1/E6E7 were grown in MEM and EMEM (MEM plus NEAA) supplemented with 10% fetal bovine serum (Gibco, NY, USA), and 1% penicillin/streptomycin (Gibco, NY, USA). Cells were grown in a humidified 5% CO₂ incubator at 37°C.

Total RNA was extracted using mirVana™ miRNA Isolation Kit (AM1561, Ambion^®). RNA sequencing (RNA-seq) was performed by Annoroad Gene Technology Co., Ltd. (Beijing, China) (Li et al., 2019). For each sample, 3 μg of total RNA were used as initial material to generate sequencing libraries using the NEBNext^® Ultra^™ Directional RNA Library Prep Kit according to the manufacturer’s recommendations. After the library was constructed, a series of processes were performed to ensure the quality of the library. Then, paired-end sequencing was performed on a single lane of Illumina HiSeq X Ten platform (Illumina, San Diego, CA, USA), with PE150 setting, producing 250 bp reads per end, according to manufacturer’s instructions.

After removing contaminated reads for adapters and low-quality reads, Bowtie2 (v2.2.6) was used for building the genome index, and clean data and was then aligned to the reference genome using HISAT2 v2.1.0. The Integrative Genomics Viewer (IGV) was used to view the mapping result (heatmap, histogram, and scatter plot). Reads count for each gene in each sample were counted by HTSeq v0.6.0, and FPKM (Fragments Per Kilobase Millon Mapped Reads) was then calculated to estimate the expression level of genes in each sample. DESeq2 v1.6.3 was designed for differential gene expression analysis between two samples with biological replicates under the theoretical basis that obeys the hypothesis of the negative binomial distribution for the value of count. DESeq2 was used to estimate the expression level of each gene in per sample by the linear regression, and then calculate the p-value with Wald test. Genes with q  ≤ 0.05 and |log₂_ratio| ≥ 1 were considered differentially expressed genes (DEGs). The Gene Ontology (GO) enrichment of DEGs was implemented by the hypergeometric test, in which the p value is calculated and adjusted as q value, and the data background is genes in the whole genome. GO terms with q  < 0.05 were considered to be significantly enriched. The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment of DEGs was implemented by the hypergeometric test, in which the p value was adjusted by multiple comparisons as q value. KEGG terms with q  < 0.05 were considered to be significantly enriched.

Hi-C sequencing (Hi-C-seq) was performed by Annoroad Gene Technology Co., Ltd. (Beijing, China). Cells were first crosslinked and then were lysed under ice condition to the extracted DNA. After quantification, the final sequencing library was diluted to 1 ng/μl. StepOnePlus™ Real-Time PCR system was used for qPCR to accurately quantify the concentration of the library. The TruSeq PE Cluster Kit v3-cbot-HS (Illumia) reagent was used to generate clusters on the cBot. After that, the library was sequenced on Illumina HiSeq X Ten platform with PE150 setting.

For mapping and interaction identification, Bowtie2 (v2.2.6) was used to map the sequenced reads to the Arabidopsis TAIR10 genome. Multiple mapped reads, unmapped paired-end reads, singleton reads, and PCR duplications were filtered by Hi-C Pro pipeline. The uniquely valid paired-end reads were kept for downstream analysis. ICE methods (Yaffe and Tanay, 2011) were used to remove different biases after building the raw contact matrices. For compartment A/B identification, the matrix2 compartment module of the cworld software was used to detect the compartment under 40 kb resolution contact matrix. Interaction distances which were below 1 MB were filtered. The lowess-smoothed average method was used to calculate the expected scores of the intra- and interaction matrix. The observed/expected ratio was log2 transformed. The patterns of chromosomal interactions were calculated at each pair of bins by using Pearson correlation, and then using this correlated matrix to do the principal component analysis. The first principal component’s eigenvalue was plotted that positive values are referred to as compartment A, which means “open chromatin”, and negative values are referred to as compartment B, also means “closed chromatin”. The gene density was defined by calculating each bin’s gene number (Lieberman-Aiden et al., 2009).

ChIP sequencing (ChIP-seq) was performed by Annoroad Gene Technology Co., Ltd. (Beijing, China). Cells were crosslinked with 1% formaldehyde for 10 min at room temperature and quenched with 125 mM glycine. The fragmented chromatin fragments were pre-cleared and then immunoprecipitated with Protein A + G Magnetic beads coupled with anti-H3K27ac (ab4729, Abcam, USA) antibody. After reverse crosslinking, ChIP and input DNA fragments were end-repaired and A-tailed using the NEBNext End Repair/dA-Tailing Module (E7442, NEB) followed by adaptor ligation with the NEBNext Ultra Ligation Module (E7445, NEB). The DNA libraries were amplified for 15 cycles and sequenced using Illumina NovaSeq 6000 with single-end 1 × 75 as the sequencing mode.

Raw reads were filtered to obtain high-quality clean reads by removing sequencing adapters, short reads (length <35 bp) and low quality reads using Cutadapt (v1.9.1) (Martin, 2011) and Trimmomatic (v0.35) (Bolger et al., 2014). Then FastQC (Kyi et al., 2018) was used to ensure high reads quality. The clean reads were mapped to the human genome (hg19) using the Bowtie2 (v2.2.6) (Langmead and Salzberg, 2012) software. Peak detection was performed using the MACS (v2.1.1) (Zhang et al., 2008) peak finding algorithm with 0.01 set as the p value cutoff. Annotation of peak sites to gene features was performed using the ChIPseeker R package (Yu et al., 2015).

The protein–protein interaction (PPI) network was constructed using Search Tool for the Retrieval of Interacting Genes (STRING) (https://www.string-db.org/) and visualized with the Cytoscape software (Shannon et al., 2003). A combined score of >0.9 was retained in the further analysis. The Gene Expression Profiling Interactive Analysis (GEPIA) was used to assess the expression of KLF12 in the three major gynecological tumors (cervical cancer, endometrial cancer, and ovarian cancer) of the TCGA data (Tang et al., 2017). The threshold Fold-change and p value were set at 1.5 and 0.05 respectively to get the expression boxplot. The Human Protein Atlas database (HPA) (https://www.proteinatlas.org) was used to analyze protein expression of KLF12 between normal and cervical cancer tissues (Uhlén et al., 2015; Tebani et al., 2020).

CTCF ChIP-seq datasets for NHEK and SiHa cells were downloaded from the NCBI Gene Expression Omnibus (GEO; accession number: GSM733740 for NHEK, and GSE143026 for SiHa) (ENCODE Project Consortium, 2012; Edgar et al., 2002).

The genes loci on the chromosomes of both SiHa and HPV16 were analyzed according to the UCSC database (http://genome.ucsc.edu/). There is a fragment of HPV16 (coordinates from 3384 to 7906/1–3132, length 7,654 bp) integrated on chromosome 13 at genomic coordinates 7378870–74087558 in the human genome ( Figure 1A ). Parts of the HPV16 E2 gene, along with the complete E4, E5, L2, L1, E6, E7, and E1 genes, were integrated into the genome of the SiHa cells. A microhomologous “AGTC” fragment was present upstream and a microhomologous “TATT” fragment was present downstream of the HPV integration. The 3’ integration locus of HPV16 is located on the second intron of LINC00393 gene, which is a high frequency integration site for HPV16 (Qiu et al., 2021), while the 5’ integration site is located on the intragenic region.

As a result, HPV16 integration leads to a recurrent pattern of DNA amplifications, with two insertional breakpoints directly flanked a twice-amplified segment ( Figure 1B ) (Akagi et al., 2014). This structural alteration is mostly represented by the rearrangements adjacent to the integration sites that leave the two HPV16 integration fragments sharing the same transcriptional orientation of LINC00393. Therefore, the analysis validated the integration status and pattern of HPV16 in SiHa cells and showed the alterations in the chromosomal structure.

To get a better understanding of the structural variations caused by HPV16 integration, we compared the Hi-C-seq data of the cervical cancer SiHa and normal human epidermal keratinocyte NHEK cells. After filtering out the same compartments (AA, BB) and the different (AB, BA) linkers, approximately 90.8 million for NHEK and 153.7 million for SiHa reads were obtained ( Table S1 ). The numbers of valid reads for NHEK and SiHa were 35,391,673 and 47,448,435, respectively. The visualized overall heatmap revealed higher order genomic organization ( Figure 2A ). The SiHa heatmap was consistent with that of the NHEK cells, except for some obvious differential organizations, which were possibly caused by cervical cancer itself. Next, we checked the A/B compartments and found that 71.2% of the total genome remained conserved between NHEK and SiHa cells, with only 13.1 and 15.7% of the annotated genome changing from A to B or B to A compartments, respectively. We also observed the A/B compartment switching at a high proportion in chromosomes such as chromosomes 2, 3, 6, and 9 ( Figure 2B ).

To investigate the differential interactions near the HPV16-LINC00393 locus, we compared the 3D structures of chromosome 13 in SiHa and NHEK cells. It showed that some genomic regions exhibited differential variations, especially around the integration locus ( Figure 2C ). Further comparison at a resolution of 1 MB disclosed several changes in the compartments composition in SiHa cells. Moreover, SiHa cells exhibited a higher degree of chromatin openness, which may facilitate gene transcriptional regulation ( Figure 2D ). A/B compartments analysis showed that 77.3% of the annotated chromosome 13 remained unchanged between the NHEK and SiHa cells, with only 8.3 and 14.4% changing from A to B or B to A compartments, respectively ( Figure 2E ).

TADs are important stable regulatory units, whose alterations are closely related to tumor development (Valton and Dekker, 2016). To identify the different TAD boundaries between SiHa and NHEK, TAD boundaries were detected at a resolution of 40 kb in each cell lines. The number (1,196 and 1,443 TAD boundaries for SiHa and NHEK, respectively) and average size of TAD boundaries were slightly different within the two cells ( Figure S1 ). Among the identified 2,332 TAD boundaries, 307 TAD boundaries were overlapped, 889 and 1,136 were unique to SiHa and NHEK, respectively ( Figure S2 ). Next, we characterized the specific changes of TAD boundaries on chromosome 13. It showed that HPV16-LINC00393 integration generated a longer chromosome, together with alterations on TAD boundaries ( Figure 2F ). After comparing the TAD boundaries around the integration site, a newly TAD boundary was observed in SiHa, along with the loss of normal interaction ( Figure 2G ). Further analysis suggested that HPV16-LINC00393 integration has induced the split of the TAD, which thereby influenced the original interaction ( Figure 2H ).

Because the TAD fusion was associated with CCCTC-binding factor (CTCF) changes, we applied ChIP-seq to detect CTCF profiles in SiHa and NHEK cells. As is shown, genome-wide CTCF distribution did not differ significantly between the two cells ( Figure 3A ). However, the binding locations of CTCF in SiHa only matched approximately 27% of the TAD boundaries in NHEK cells ( Figure 3C ). As for the HPV-LINC00393 integrated chromosome 13, the number of CTCF was 3,214 and 5,635, respectively, with 1,322 CTCF being shared between the two cells ( Figure 3E ). The number of CTCF targeted genes located in the seven different gene regions was also discrepant ( Figure 3G and Table S2 ). Furthermore, we confirmed that SiHa cells displayed intense CTCF binding at the TAD boundary of LINC00393 locus, with a weak binding found in the NHEK cells ( Figure 3H ).

Since histone epigenetic modifications can also disturb host 3D genome structure, we then performed H3K27ac ChIP-seq in the two cells. H3K27ac distribution was narrow in SiHa cells, indicating reduced enhancer activity compared to NHEK cells ( Figure 3B ). Meanwhile, the binding locations of H3K27ac in SiHa matched approximately 40% of that in NHEK ( Figure 3D ). Using a q value cutoff of 0.05, we found that approximately 48% of the H3K27ac regions in SiHa and approximately 29% of the H3K27ac peaks in NHEK were overlapped ( Figure 3F ). The number of H3K27ac modified genes was decreased in the seven different gene regions of SiHa cells ( Figure 3G and Table S2 ). Furthermore, it showed that SiHa cells exhibited weak H3K27ac binding at the TAD boundary of the LINC00393 integration site, whereas an intense binding was observed in the NHEK cells ( Figure 3I ).

Chromatin translocation was reported to alter host gene expression and lead to cancer development (Gryder et al., 2021). To investigate the effect of translocation induced changes of TAD boundary on the expression of surrounding genes, we performed the transcriptome analysis of SiHa and normal cervical cell Ect1/E6E7. The overall distribution of the differentially expressed genes (DEGs) was shown in Figures 4A–C . Next, we investigated the effect of LINC00393 integration on the expression of genes present on chromosome 13. A total of 74 genes were characterized as DEGs, with 37 downregulated and 37 upregulated ( Figure 4D andd Table S3 ). The nearby regions of LINC00393 such as KLF5, KLF12, TBC1D4, MYCBP2, and SCEL showed downregulated expression, while COMMD6 appeared as an upregulated gene ( Figure 4E ). Further analysis revealed that these DEGs were potentially involved in tumorigenesis-related, DNA repair, and HPV infection pathways ( Figure 4F ). Because transcription factor (TF) is an important determinant responsible for transcriptional activation of genes involved in tumorigenesis, we conducted TF analysis based on the DEGs. As is shown, the most enriched TF belonged to C2H2 type zinc finger (zf-C2H2) family ( Figure 4G ).

The tumor suppressor gene KLF12 is a member of zf-C2H2. However, whether KLF12 plays a direct role in cervical cancer and whether its inhibition can promote cervical tumorigenesis remains unclear. To address this, we used the STRING database and mapped the protein–protein interaction network of KLF12 ( Figure 5A ). Among the top ten proteins that interact with KLF12, CTBP1 was reported to be a key mediator for the transcriptional inhibitory role of KLF12 (McConnell and Yang, 2010). We next used the GEPIA website based on the TCGA database and analyzed the expression of KLF12 in the three major gynecological tumors. Compared with the healthy women, KLF12 expression was significantly downregulated in patients with ovarian cancer, endometrial cancer, and cervical cancer ( Figure 5B ). Furthermore, we used the HPA database to validate KLF12 expression in cervical cancer and normal cervical tissues. Decreased expression of KLF12 was observed in the nucleus of both cervical squamous cell carcinoma tissue and adenocarcinoma tissue ( Figure 5C ). These data indicate that KLF12 may connect with cervical tumorigenesis and have the potential to predict the onset of cervical cancer in HPV16-LINC00393 integrated patients.