Comprehensive Analysis of circRNA Expression Profiles During Cervical Carcinogenesis

Circular RNAs (circRNAs) are regulatory molecules that participate in the occurrence, development and progression of tumors. To obtain a complete blueprint of cervical carcinogenesis, we analyzed the temporal transcriptomic landscapes of mRNAs and circRNAs. Microarrays were performed to identify the circRNA and mRNA expression profiles of cervical squamous cell carcinoma (CSCC) and high-grade squamous intraepithelial lesion (HSIL) patients compared with normal controls (NC). Short time-series expression miner (STEM) was utilized to characterize the time-course expression patterns of circRNAs and mRNAs from NC to HSIL and CSCC. A total of 3 circRNA profiles and 3 mRNA profiles with continuous upregulated patterns were identified and selected for further analysis. Furthermore, functional annotation showed that the mRNAs were associated with DNA repair and cell division. The protein-protein interaction (PPI) network analysis revealed that the ten highest-degree genes were considered to be hub genes. Subsequently, a competing endogenous RNA (ceRNA) network analysis and real-time PCR validation indicated that 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 were highly correlated with cervical carcinogenesis. Silencing of hsa_circ_0003954 inhibited SiHa cell proliferation and perturb the cell cycle in vitro. This study provides insight into the molecular events regulating cervical carcinogenesis, identifies functional circRNAs in CSCC, and improves the understanding of the pathogenesis and molecular biomarkers of CSCC and HSIL.


INTRODUCTION
Cervical cancer is one of the most prevalent gynecological malignancies and the fourth most fatal cancer (1). Cervical squamous cell carcinoma (CSCC) is the most frequent pathological subtype of cervical cancer, accounting for approximately 80-85% of cases (2). Uterine cervix carcinogenesis is a step-by-step process that shows a continuum of neoplastic transitions from the persistence of high-risk human papillomavirus (HR-HPV, mostly HPV16) infection to low-grade squamous intraepithelial lesion (LSIL) to high-grade squamous intraepithelial lesion (HSIL) and then to invasive cancer histologically (3). It has been reported that almost half of HSILs will progress to invasive cancer (4,5). However, the molecular mechanism of cervical carcinogenesis is not completely understood, and previous studies of cervical precancerosis have been extremely limited. Hence, exploration of the genetic changes and epigenetic modifications may reveal new clues to elucidate the molecular mechanisms of cervical carcinogenesis.
Circular RNAs (circRNAs), a kind of noncoding RNA without a 3′ tail or 5′ cap, are naturally occurring endogenous molecules (6,7). The microRNA (miRNA) response elements (MREs) in circRNAs can suppress the function of miRNAs, during this process, circRNAs inhibit the activity of miRNAs and modulate the expression of their downstream target genes in various types of malignancies (8). Previous research has highlighted the important roles of circRNAs in cancer occurrence and progression (9)(10)(11)(12)(13)(14).
Recently, multiple studies have identified the participation of circRNAs in CSCC formation and progression (15,16). Cervical cancer tissue samples and cell lines both highly expressed CiRS-7, and its expression levels correlated significantly with prognosis (17). CircE7 in CaSki cervical carcinoma cells inhibits the growth of cancer cells by reducing the protein levels of E7 (18). However, while the majority of research has primarily focused on identifying the differentially expressed circRNAs (DECs) between cervical cancer and controls, few researchers have focused on HSIL. The role of DECs in cervical carcinogenesis has not been systematically described.
In this study, we recruited not only CSCC patients and controls but also HSIL individuals. Due to the stepwise progression of CSCC, the expression of circRNAs and mRNAs could be dysregulated at any specific tumorigenesis stage. To characterize the changes in circRNA and mRNA expression, we performed trend analysis to identify the predominant circRNAs and mRNAs in the control, HSIL and CSCC groups. This study provides a starting point for further research into the molecular mechanism of circRNAs in cervical carcinogenesis, which provides new insight into the multiple and complex factors in early cervical carcinogenesis.

Patients and Specimens
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 HPV16positive 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.

RNA Extraction and Quality Control
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 and Data Analysis
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).

Functional Enrichment Analysis
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.

Short Time-Series Expression Miner (STEM)
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).

Protein-Protein Interaction (PPI) Network Construction
The PPI network of DEMs was constructed using the STRING (a search tool for the retrieval of interacting genes, http://stringdb.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.

qRT-PCR and CSCD (Cancer-Specific CircRNA Database) Analysis
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 -DDCt 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.

Transfection
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 and Cell Cycle Analysis
Cell proliferation was detected through the CCK-8 assay (Meilunbio, Dalian, China). For transient transfection experiments, 1×10 3 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.

Statistical Analysis
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).

Analysis of circRNA and mRNA Expression Profiles
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.

Functional Annotation of DEMs
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, ubiquitinmediated proteolysis, protein processing in endoplasmic reticulum, p53 signaling pathway and cell cycle ( Figure 3D).

Temporal Gene Expression Patterns of mRNAs and circRNAs
To explore the differences in gene expression between cervical carcinogenesis stages, the STEM tool was applied to profile stagespecific 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.

The PPI Network
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.

The ceRNA Network
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.

Silencing of hsa_circ_0003954 Effects Cell Proliferation Progression and Cell Cycle in SiHa Cells
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 The blue circles represent protein-coding mRNAs, the yellow triangles represent miRNAs and green rhombuses represent circRNAs. The gray solid lines represent the circRNA-miRNA regulatory relationships. The size of the node represents its degree in the network, and the larger the node is, the higher the degree. ceRNA, competing endogenous RNA. 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).

DISCUSSION
CircRNAs, playing in a variety of human diseases, contributes to the pathogenesis of abundant cancers due to dysregulated expression, including CSCC (27)(28)(29)(30)(31)(32)(33)(34). Nevertheless, because the dynamics of gene expression are characterized by a phasic pattern and cervical carcinogenesis is a gradual process (35), the role of circRNAs in the development of CSCC cannot be fully characterized.
In the present study, we provided comprehensive profiling of the transcriptome involving circRNA and mRNA from NC to HSIL and CSCC. What's more, we identified 3 mRNA profiles and 3 circRNA profiles that were expressed in a time-dependent manner during the carcinogenic process of the cervix. Interestingly, we observed a significant enrichment of cancerrelated signaling pathways in the DEMs over time, such as DNA replication, cell cycle, and DNA repair-related GO functions, based on the DEMs obtained by STEM. Followed by a PPI network ( Figure 6), CDK1, BUB1, KIF11, NDC80, BUB1B, CCNB2, PCNA, CCNB1, MAD2L1 and CDCA8 were found to be hub genes with an extremely high degree of intramodular connectivity. Among these hub genes, CDK1, KIF11, CCNB2, PCNA and CCNB1 were found in the two hub modules identified by MCODE. Previous studies have identified three genes, CDK1, PCNA and CCNB1 that contribute to the occurrence and development of CSCC (36,37). However, there is no evidence linking BUB1 and KIF11 with CSCC or of their association with circRNAs. Overall, these hub genes are relevant to carcinogenesis, and further study of their role in CSCC is warranted.
Further alignments and prediction of the selected DECs and DEMs led to the establishment of a ceRNA regulatory network ( Figure 7). Notably, 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. It has been found that hsa_circ_0001955/miR-145-5p is the key axis in the carcinogenesis of colorectal cancer (38) and in vitro, hsa_circ_0001955 promotes hepatocellular carcinoma cell proliferation, migration, and invasion via miR-145-5p/NRAS (39). However, little is known about the role of hsa_circ_0001955 and hsa_circ_0003954 in CSCC. CDK1, target genes of hsa_circ_0001955, have been reported to participate in cervical cancer proliferation, invasion and migration and have been associated with tumor stage and lymph node status (40). Additionally, SYCP2, another important gene identified in our analysis, may contribute to HPV-associated cancer development, according to recent research (41). However, no studies have found that NEDD4L plays an essential regulatory role in CSCC. It has been reported that miR-15a-3p increased radiosensitivity in cervical cancer by targeting TPD52, suggesting that miR-15a-3p may be a potential therapeutic target (42). Little is known about how miR-6719-3p and hsa-miR-1277-5p influence the development and progression of cervical cancer.
To determine the roles of circRNA in the carcinogenesis of cervix, molecular biological experiments are necessary, since our results were based on computational prediction. According to the biological functions and potential mechanisms in cervical carcinogenesis, DNA replication and cell cycle play important roles in CSCC. Therefore, we performed related experiments and found that knocking down hsa_circ_0003954, the top upregulated circRNA, inhibits SiHa cell proliferation and modulates the cell cycle in vitro.
In summary, we characterized the circRNA and mRNA transcriptomes of cervical tissues during cervical carcinogenesis. We used STEM analysis to determine the temporal patterns of circRNA The results are presented as the mean ± SD (NC-qRT-PCR, n = 5; HSIL-qRT-PCR, n = 5; CSCC-qRT-PCR, n = 5). **p < 0.01, ***p < 0.001, ****p < 0.0001. SD, standard deviation. and mRNA in cervical carcinogenesis, and we also identified two extremely highly expressed circRNAs, hsa_circ_0001955 and hsa_circ_0003954, that play a potential role in cervical carcinogenesis. Overall, the results of our study shed light on the molecular mechanisms, providing new evidence and insights of CSCC.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi. nlm.nih.gov/, GSE166466. (E-G) The cell cycle progression was analyzed by flow cytometry after transfected with si-circ or si-NC. ***p < 0.001, ****p < 0.0001.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Ethics Committee of Second Hospital of Shanxi Medical University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
WW designed and supervised this project. HL, YL, and WW analyzed the data, and wrote the manuscript. HL and YL performed the experiments, YZ, JC, XZ, BZ, LG, and WW revised the manuscript. HL and YL contributed to data interpretation. All authors contributed to the article and approved the submitted version.