A total of 35 human cervical specimens were acquired, including 12 patients with only HPV16-positive normal controls (NC), 11 with only HPV16-positive HSIL and 12 with only HPV16-positive stage IA-IIA CSCC, who were collected at the Second Hospital of Shanxi Medical University between December 2019 and October 2020. We collected all patients’ clinical information, including age, HPV testing, the ThinPrep^® Pap Test (TCT) and pathology (by colposcopy biopsy, cervical conization or hysterectomy) results. Based on the International Federation of Gynecology and Obstetrics (FIGO) criteria, the clinical stage of the patients with CSCC was determined (19). The diagnosis of all of the cases was histologically confirmed by two independent pathologists, and all of the CSCC ( Figures S1A–C ) and HSIL ( Figures S1D–F ) samples were assessed by hematoxylin and eosin (H&E) staining. Cases were excluded when they met any of the following exclusion criteria: 1) had a history of cervical physical therapy (ablation and cryosurgery), chemotherapy and pelvic radiotherapy; 2) had other tumors; and 3) had an immunocompromised condition (e.g., infection with human immunodeficiency virus). The clinicopathological characteristics of the patients and the usage of specimens are summarized in Table S1 . This study was approved by the Ethics Committees of the Second Hospital of Shanxi Medical University.

SiHa (human cervical squamous cell carcinoma cell line, HPV 16-positive) and HcerEpic (human normal cervical epithelial cell line) cells (ATCC, Rockville, MD, USA) were cultured with DMEM (General Electric Company, USA) with 10% fetal bovine serum (Gibco; Thermo Fisher Scientific, Inc., Waltham, MA, USA), 100 mg/mL streptomycin, and 100 U/mL penicillin in a humidified atmosphere of 5% CO₂ and 37°C.

Total RNA was extracted from tissues using TRIzol Reagent (Life Technologies, CA, US) and purified with an RNeasy Mini Kit (Qiagen, Valencia, CA, USA). Assaying the purity and integrity of the RNA was achieved using agarose gel electrophoresis and a UV/vis spectrophotometer (Thermo, NanoDrop 2000, USA).

Microarray hybridization was performed in accordance with the Affymetrix GeneChip Expression Analysis Technical Manual (www.affymetrix.com). The data were analyzed using the robust multichip analysis (RMA) algorithm using default Affymetrix settings. The values presented are the log2 RMA signal intensity. Differentially expressed mRNAs (DEMs) and DECs among the three groups were filtered according to the following criteria: p < 0.05 and fold change > 1.2. The microarray data were uploaded to a public database (accession number: GES166466).

The online tool Metascape (http://metascape.org/) was used to analyze DEMs (20). Three types of gene ontology (GO) were used to functionally enhance biological process (BP), molecular function (MF) and cellular component (CC) (21). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment was also performed.

STEM software version 1.3.13 (http://www.cs.cmu.edu/~jernst/stem) was utilized to cluster and view probable circRNA and mRNA expression patterns over time (22). The STEM clustering method and other options were set as defaults. The circRNA and mRNA expression profiles were clustered based on statistically significant values (p < 0.05).

The PPI network of DEMs was constructed using the STRING (a search tool for the retrieval of interacting genes, http://string-db.org) online database (23). In the current study, a combined score ≥ 0.99 was considered the cutoff criterion. The PPI network was visualized on a free bioinformatics platform provided by Cytoscape version 3.7.1 (https://cytoscape.org/). MCODE (molecular complex detection) version 3.7.1 The Cytoscape plugin was used to screen the potential hub modules.

CircRNA–miRNA interactions were predicted using the miRanda (http://www.microrna.org/microrna/home.do) online database (24), and miRNA–mRNA interactions were predicted using the TargetScan (http://www.targetscan.org/vert_72/) (25) and miRanda prediction databases. The ceRNA networks were visualized by Cytoscape.

Quantitative real-time polymerase chain reaction (qRT-PCR) was performed using Q SYBR Green Supermix (Bio-Rad, Hercules, CA, USA), and PCR-specific amplification was conducted in the 7900 HT Sequence Detection System (ABI PRISM; Waltham, MA, USA). The expression was determined by using the threshold cycle (Ct) method, and relative expression levels were calculated via the 2^–ΔΔCt method. The top five circRNAs with the highest degree in the ceRNA network and the top five miRNAs and mRNAs for the potential circRNAs were validated by qRT-PCR. The primers for these RNAs are presented in Table S2 . CSCD (http://gb.whu.edu.cn/CSCD/) (26), an online tool to investigate cancer-specific circRNAs, was utilized to obtain the structure of potential circRNAs.

To transfect with siRNA, we used custom-designed siRNAs targeting hsa_circ_0003954 ( Table S3 and Figure 10B ). For transfection, SiHa cells were grown on 6 well plates. Cells were transfected at 24 h with 30 pmol siRNA or scrambled control (GenePharma, Shanghai, China) using Lipofectamine 3000 (Invitrogen MA, USA) according to the manufacturer’s protocol. A total of three biological triplicates were conducted.

Cell proliferation was detected through the CCK-8 assay (Meilunbio, Dalian, China). For transient transfection experiments, 1×10³ cells were plated in 96 well plates for 24 hours at 37°C. Proliferation absorbance was measured with a multifunctional microplate reader (SpectraMax M5, MD, USA). Experiments were repeated three times. Cell cycle assays were conducted using propidium iodide stained SiHa cells by a Beckman Coulter FC500 flow cytometer (Beckman-Coulter, Hialeah, FL) and analyzed using Modfit software.

Experimental data are presented as the mean ± standard deviation (SD) of at least three experiments. Significant differences were assessed by Student’s t-test. P < 0.05 was considered statistically significant. Figures were drawn using R Studio version 3.3.4 (The R Foundation for statistical computing, Vienna, Austria).

Genome-wide analysis of differentially expressed profiles in circRNA and mRNA among CSCC, HSIL and NC tissues was performed. A total of 3172 circRNA candidates were differentially expressed with a fold change ≥ 1.2 (p < 0.05). Thirty-two circRNAs revealed a fold change ≥ 10. Hsa_circ_0066984 (fold change ~ 24) was the most dysregulated circRNA. The candidate DECs were distributed on 46 human chromosomes, including the chromosomes 1, 2, 3 and X chromosome, which contained more circRNAs than the other chromosomes ( Figure 1 ). A total of 4519 mRNAs among CSCC, HSIL and NC tissues had a fold change ≥ 1.2 (p < 0.05). Summarization of the coding gene expression profile showed that 46 mRNAs displayed a fold change ≥ 10. CircRNA and mRNA expression patterns among CSCC, HSIL and NC samples were significantly differentially expressed, as shown by hierarchical clustering ( Figure 2 ). The cluster heatmaps of the DECs and DEMs showed good discrimination among CSCC, HSIL and NC samples.

GO and KEGG pathway enrichment analyses were conducted to explore potential biological processes and pathways enriched by DEMs. Figure 3 presents the top ten enriched BP, CC, MF terms and KEGG pathways. The enriched BP terms were mainly related to carcinogenic processes, such as viral transcription, translational initiation, regulation of mRNA stability and positive regulation of ubiquitin−protein ligase activity involved in regulation of mitotic cell cycle transition ( Figure 3A ), the most enriched CC was the nucleus ( Figure 3B ) and the most enriched MF was protein binding ( Figure 3C ). Similar to the GO term results, KEGG pathway enrichment analysis found that DEMs were mainly enriched in viral carcinogenesis, ubiquitin-mediated proteolysis, protein processing in endoplasmic reticulum, p53 signaling pathway and cell cycle ( Figure 3D ).

To explore the differences in gene expression between cervical carcinogenesis stages, the STEM tool was applied to profile stage-specific gene expression patterns. The microarray data were normalized to the NC data, and the temporal gene expression profiles were identified. Six temporal mRNA profiles and six temporal circRNA profiles were statistically significant (p < 0.05), including profiles 10, 11, 12, 13, 14 and 15 ( Figure 4 ). Profiles 12, 13 and 15 had the same continuous upregulation patterns ( Figures 5A–C ) and were selected for further analysis, which mostly related to the cell cycle ( Figures 5D, E ), DNA replication ( Figure 5D ), DNA repair ( Figure 5F ) and cell division ( Figures 5D–F ) according to the BP analysis.

After removing the unconnected nodes and nodes that could not connect to the main network, a PPI network consisting of 267 interactions among 95 nodes was constructed ( Figure 6A ). The top ten genes with the highest degree of connectivity, namely, CDK1, BUB1, KIF11, NDC80, BUB1B, CCNB2, PCNA, CCNB1, MAD2L1 and CDCA8, were selected as hub genes. MCODE in Cytoscape was applied to identify hub modules in the PPI network, which revealed the biological functions of the key protein complexes with the highest degree of connectivity in HSIL and CSCC tissues. In addition, Figures 6B, C illustrate that the top two modules with the highest score were selected as potential hub modules (Modules 1 and 2), where hub genes such as CDK1, KIF11, CCNB2, PCNA and CCNB1 were included.

A total of 184 target miRNAs for DECs and DEMs were obtained from the miRanda and TargetScan databases. We constructed a ceRNA network to further investigate the role of mRNAs and circRNAs in cervical carcinogenesis. As shown in Figure 7 , a total of 479 interactions between the selected genes were identified and visualized. Multiple circRNAs could act as ceRNAs to capture downstream miRNAs, thus influencing the phenotype by regulating mRNAs. Furthermore, hsa-miR-1277-5p, hsa-miR-335-3p, hsa-miR-153-5p, hsa-miR-30a-3p and hsa-miR-412-3p were identified as the miRNAs with the highest degree in the network.

qRT-PCR analysis of CSCC, HSIL and NC samples was performed to validate the top 5 circRNAs, namely, hsa_circ_0016456 (degree=64), hsa_circ_0008617 (degree=30), hsa_circ_0001955 (degree=26), hsa_circ_0003954 (degree=21) and hsa_circ_0076726 (degree=14) ( Table 1 ). The expression patterns of hsa_circ_0001955 and hsa_circ_0003954 were both continuously upregulated in NC to HSIL and CSCC by qRT-PCR (p < 0.01), which was consistent with the above microarray analysis results ( Figure 8 ). Moreover, the top 5 miRNAs and mRNAs predicted by ceRNA analysis for hsa_circ_0001955 and hsa_circ_0003954 ( Figures 9C, D ) were also validated by qRT-PCR ( Figures 9E, F ). The expression patterns of hsa_circ_0001955/hsa-miR-6719-3p/CDK1, hsa_circ_0001955/hsa-miR-1277-5p/NEDD4L and hsa_circ_0003954/hsa-miR-15a-3p/SYCP2 fit with ceRNA network by qRT-PCR (p < 0.01). The structures of hsa_circ_0001955 and hsa_circ_0003954 are presented ( Figures 9A, B ) based on the data from CSCD, which indicated that both circRNAs contained MREs.

Hsa_circ_0003954, the top upregulated circRNA validated by qRT-PCR, was selected as the prospective circRNA for further research. The expression of hsa_circ_0003954 was evaluated by qRT-PCR in SiHa and HcerEpic cell lines, and the results showed that hsa_circ_0003954 expression was higher in SiHa cells than in HcerEpic cells ( Figure 10A ). Three siRNAs targeting the junction sites of hsa_circ_0003954 ( Figure 10B ) were designed and the qRT-PCR results showed that the expression of hsa_circ_0003954 was significantly downregulated in SiHa cells transfected with the siRNA segments ( Figure 10C ). Of the three siRNAs, si-circRNA#2 was selected for further investigation with the highest silencing effectiveness in SiHa cells.

In the proliferation assay ( Figure 10D ), growth curves performed by CCK-8 assays demonstrated that silencing hsa_circ_0003954 significantly inhibited the proliferation viability of SiHa cells (p < 0.001). The cell cycle analysis showed that more SiHa cells were distributed in G1 phase and less in S phase after silencing hsa_circ_0003954, which suggested that SiHa cells were arrested at G1 phase by silencing hsa_circ_0003954 ( Figures 10E–G ).