Abstract Purpose Cervical squamous cell carcinoma (CSCC) seriously affects women’s health worldwide, and it is of great significance to illuminate the specific role of circRNAs in CSCC. Materials and Methods Three mRNA datasets, two miRNA datasets and one circRNA dataset of CSCC, downloaded from GEO, were utilized in this study. Differentially expressed mRNAs (DEmRNAs), miRNAs (DEmiRNAs) and circRNAs (DEcircRNAs) were identified, and a ceRNA (DEcircRNA-DEmiRNA-DEmRNA) regulatory network was constructed. GO and pathway analyses of DEcircRNAs and DEmRNAs in the ceRNA regulatory network were performed. Quantitative real-time polymerase chain reaction (qRT-PCR) validation of the expression of the selected DEmRNAs, DEmiRNAs and DEcircRNAs was performed. Results A total of 1356 DEmRNAs, 13 DEmiRNAs and 77 DEcircRNAs were obtained. The ceRNA network contained 3 circRNA-miRNA pairs and 158 miRNA-mRNA pairs, including 3 circRNAs, 3 miRNAs, and 138 mRNAs. Functional annotation of DEmRNAs in the ceRNA regulatory network revealed that these DEmRNAs were significantly enriched in cell cycle, p53 signalling pathway and DNA replication. The qRT-PCR results were generally consistent with those of our integrated analysis. Conclusion In conclusion, we speculate that the regulation of the hsa_circ_0000069/hsa-miR-125b-5p/CDKN2A and hsa_circ_0020594/hsa-let-7c-5p/CCNB2 axes may be involved in CSCC. Keywords: cervical squamous cell carcinoma, CSCC, circRNA, ceRNA network, GEO Introduction Cervical cancer (CC) is one of the most widespread gynaecological cancers that threatens women’s health in developing countries.[36]^1 CC progresses through a long stage of precancerous lesions, called cervical intraepithelial neoplasia (CIN), which is classified into three grades: CIN 1, CIN 2 and CIN 3.[37]^2 Human papillomaviruses (HPVs) have been widely considered the main cause of CC.[38]^3 Cervical squamous cell carcinoma (CSCC) is the most frequent type, accounting for the majority of CCs.[39]^4 The molecular mechanisms of CSCC are complex and have not yet been fully elucidated. Hence, it is necessary to deeply understand the molecular mechanisms and pathways of CSCC progression. In contrast to traditional linear RNAs, circRNAs are covalently closed circular molecules that have attracted wide attention from researchers.[40]^5 CircRNA, as one type of competing endogenous RNA (ceRNA), can bind to miRNAs, acting as a “miRNA sponge”, and regulate their target genes at the transcriptional or post-transcriptional level, thereby affecting the biological behaviour of tumours.[41]^6 The advent of microarray-based technology has facilitated studies on the molecular mechanisms of different tumours.[42]^7 However, investigations on CSCC are still limited, particularly with regard to circRNAs in CSCC. Wang et al suggested that hsa_circ_0101996 combined with hsa_circ_0101119 can serve as potential biomarkers for CSCC detection.[43]^8 Yi et al identified five circRNAs that may play critical roles in CC.[44]^1 Ma et al revealed that circRNA-000284 is involved in cell proliferation and invasion of CC by sponging miR-506.[45]^9 Therefore, more CSCC-related circRNAs need to be unearthed. The current study performed an integrated analysis with datasets of mRNA, miRNA and circRNA expression profiles of CSCC collected from the GEO database. The identification of differentially expressed mRNAs (DEmRNAs), miRNAs (DEmiRNAs) and circRNAs (DEcircRNAs) and the construction of a ceRNA (DEcircRNA-DEmiRNA-DEmRNA) regulatory network were conducted. This present work sought to help elucidate the underlying mechanisms of CSCC pathogenesis and provide additional knowledge for CSCC. Materials and Methods Dataset Collection To acquire the mRNA, miRNA and circRNA expression profiles of CSCC tumour tissues and normal cervix tissues, datasets in the GEO database with the following criteria were retrieved: selected datasets should be whole-genome mRNA/miRNA/circRNA expression profiles by array; these data were derived from CSCC tumour tissues and normal cervix tissues; (3) datasets were normalized or original. Three mRNA datasets, two miRNA datasets and one circRNA dataset were enrolled in this study ([46]Table 1). Table 1. List of mRNA/miRNA/circRNA Study Samples from GEO GEO ID Author Platform Samples (N: CSCC) Year mRNA  [47]GSE63514 den Boon J [48]GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array 24: 28 2015  [49]GSE9750 Vundavalli M [50]GPL96 [HG-U133A] Affymetrix Human Genome U133A Array 24: 33 2008  [51]GSE7803 Kuick R [52]GPL96 [HG-U133A] Affymetrix Human Genome U133A Array 10: 21 2007 miRNA  [53]GSE30656 Sie D [54]GPL6955 Agilent-016436 Human miRNA Microarray 1.0 (Feature Number version) 10: 10 2012  [55]GSE19611 Pereira PM [56]GPL7534 National DNA microarray facility of University of Aveiro miRNA chip v1.1 23: 4 2010 circRNA  [57]GSE102686 Jiao J [58]GPL19978 Agilent-069978 Arraystar Human CircRNA microarray V1 5: 5 2017 [59]Open in a new tab Differential Expression Analysis The raw data were visualized with the hist function, normalized with the log function, and centralized and standardized with the scale function. We calculated differential expression p-values and effect sizes from data either from classical or moderated t-tests (Limma, SMVar) for study and combined these p-values by the inverse normal method. Multiple comparison correction was performed by the Benjamini and Hochberg method to obtain the false discovery rate (FDR). By default, the FDR is controlled at 5%. MetaMA was utilized to acquire the DEmRNAs and DEmiRNAs with FDR < 0.05 and |Combined.ES| > 1. The Limma package in R was applied to acquire the DEcircRNAs with FDR < 0.05 and |log2 fold change| > 1. With the R package “pheatmap”, hierarchical clustering analysis of DEmRNAs, DEmiRNAs and DEcircRNAs was performed. GO and Pathway Analysis of Host Genes of DEcircRNAs GeneCodis3 ([60]http://genecodis.cnb.csic.es/analysis) was applied to perform GO and KEGG pathway enrichment analysis of host genes of DEcircRNAs. The threshold was set as FDR < 0.05. Construction of the ceRNA (DEcircRNA-DEmiRNA-DEmRNA) Regulatory Network According to the results of the differential expression analysis, DEcircRNA–DEmiRNA interaction pairs were calculated with starBase v3.0 ([61]http://starbase.sysu.edu.cn). The targeted DEmRNAs of DEmiRNAs were predicted with six bioinformatic algorithms (RNA22, miRanda, miRDB, miRWalk, PICTAR2 and TargetScan). Then, with miRWalk, the confirmed targeted mRNAs of miRNAs were obtained. Finally, the confirmed DEmiRNA-DEmRNA pairs derived from miRWalk and the DEmiRNA-DEmRNA pairs predicted by ≥4 algorithms were retained for network construction. The ceRNA (DEcircRNA-DEmiRNA-DEmRNA) regulatory network was constructed by combining circRNA-miRNA pairs with miRNA-mRNA pairs. Ultimately, Cytoscape 3.5.0 ([62]http://cytoscape.org/) was used to visualize the regulatory network. With GeneCodis3, GO and KEGG pathway analyses of all DEmRNAs in the ceRNA regulatory network were performed. Statistical significance was defined as FDR < 0.05. Validation in Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR) and TCGA Database Six paired CSCC tumour tissues and corresponding adjacent non-tumour tissues were obtained from 6 patients with CSCC. The clinical characteristics of the individuals included in this study are displayed in [63]Table S1. We obtained written informed consent from every participant. This study was approved by the ethics committee of Changning Maternity and Infant Health Hospital and in accordance with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. Total RNA was isolated with the TRIzol reagent (Invitrogen, USA) following the manufacturer’s protocol. The qRT-PCRs were performed in an ABI 7300 Real-time PCR Detection System with SuperReal PreMix Plus (Invitrogen, USA). Relative gene expression was analysed by the 2-ΔΔCT method. Human GAPDH and ACTB were used as endogenous controls for mRNA expression in the analysis. Human U6 and GAPDH were used as endogenous controls for miRNA and circRNA expression in the analysis, respectively. The qRT-PCR was performed in triplicate. The p-value was assessed by one-way analysis of variance (ANOVA), and p < 0.05 was considered the criterion of statistical significance. The mRNA and miRNA profiles of 252 CSCC cases and 2 adjacent non-tumour controls derived from the TCGA database were downloaded to validate the expression levels of selected DEmRNAs and DEmiRNAs. The difference in expression levels is displayed in box plots. In addition, to further investigate the prognostic value of selected DEmiRNAs and DEmRNAs, survival analysis was performed by using the survival package in R. Results DEmRNAs, DEmiRNAs and DEcircRNAs in CSCC Tumour Tissues Compared to Normal Cervix Tissues A total of 1356 DEmRNAs (750 up- and 606 down-regulated mRNAs) and 13 DEmiRNAs (4 up- and 9 down-regulated miRNAs) were identified in CSCC compared to normal cervix tissues with FDR < 0.05 and |Combined.ES| > 1 ([64]Figure 1A and [65]B). Among them, CDKN2A and UPK1A were the most up-regulated and down-regulated mRNAs; hsa-miR-106b-5p and hsa-miR-125b-5p were the most up-regulated and down-regulated miRNAs, respectively ([66]Table 2). A total of 77 DEcircRNAs (45 up- and 32 down-regulated circRNAs) were obtained in CSCC with FDR < 0.05, |log2 fold change| > 1 ([67]Figure 1C). Among them, hsa_circRNA_104052 and hsa_circRNA_103384 were the most up-regulated and down-regulated circRNAs, respectively ([68]Table 3). Figure 1. [69]Figure 1 [70]Open in a new tab The heatmap of top 100 up- and down-regulated DEmRNAs (A), DEmiRNAs (B) and DEcircRNAs (C) between CSCC and normal cervix tissues. Table 2. Top 10 Up- and Down-Regulated DEmRNAs and DEmiRNAs in CSCC mRNA Combined.ES p-value FDR Regulation CDKN2A 4.09523096 0 0 Up KNTC1 3.415615858 0 0 Up MCM2 3.262138669 0 0 Up DTL 3.12448687 0 0 Up ECT2 2.932967572 0 0 Up SPAG5 2.87189987 0 0 Up RFC4 2.86851678 0 0 Up MCM6 2.863540309 0 0 Up KIF14 2.834991169 0 0 Up TIMELESS 2.811416252 0 0 Up UPK1A −3.938910321 0 0 Down CRNN −3.484057128 0 0 Down ENDOU −3.480046127 0 0 Down CRISP3 −3.397031566 0 0 Down EDN3 −3.051176288 0 0 Down MAL −2.902463244 0 0 Down BBOX1 −2.892847988 0 0 Down SLURP1 −2.824556627 0 0 Down PPP1R3C −2.691069042 0 0 Down DSG1 −2.626794273 0 0 Down miRNA  hsa-miR-106b-5p 1.70322226 5.51E-05 0.00282703 Up  hsa-miR-21-5p 1.314433323 0.000352172 0.009039083 Up  hsa-miR-18a-5p 1.029429377 0.003619066 0.047549331 Up  hsa-miR-142-3p 1.019333791 0.003999802 0.047549331 Up  hsa-miR-125b-5p −1.753027386 6.34E-06 0.000507467 Down  hsa-miR-203a-3p −1.73909463 6.59E-06 0.000507467 Down  hsa-miR-99a-5p −1.414419027 0.00011931 0.004593453 Down  hsa-miR-497-5p −1.387732202 0.000150657 0.004640232 Down  hsa-let-7c-5p −1.168987104 0.001015067 0.022331465 Down  hsa-miR-195-5p −1.129936328 0.001989935 0.038306241 Down  hsa-miR-100-5p −1.059650819 0.002817986 0.047549331 Down  hsa-miR-149-5p −1.040611481 0.003198146 0.047549331 Down  hsa-miR-199b-5p −1.010825035 0.004013905 0.047549331 Down [71]Open in a new tab Abbreviations: DEmRNAs, differentially expressed mRNAs; DEmiRNAs, differentially expressed miRNAs; ES, effect size; FDR, false discovery rate; CSCC, cervical squamous cell cancer. Table 3. Top 10 Up- and Down-Regulated DEcircRNAs in CSCC circRNA Alias Log[2]FC p-value FDR Regulation hsa_circRNA_104052 hsa_circ_0008285 1.72117645 3.74E-07 0.000374916 Up hsa_circRNA_001846 hsa_circ_0000520 1.655653091 3.29E-07 0.000374916 Up hsa_circRNA_000543 hsa_circ_0000326 1.429575197 1.18E-06 0.000707729 Up hsa_circRNA_104852 hsa_circ_0006174 1.533928394 3.78E-06 0.001200045 Up hsa_circRNA_103803 hsa_circ_0072008 1.437754354 3.79E-06 0.001200045 Up hsa_circRNA_103450 hsa_circ_0002768 1.331734139 2.82E-06 0.001200045 Up hsa_circRNA_000585 hsa_circ_0000515 1.27949496 4.39E-06 0.001200045 Up hsa_circRNA_100146 hsa_circ_0011385 1.161235276 3.79E-06 0.001200045 Up hsa_circRNA_103948 hsa_circ_0003528 1.010710202 1.88E-05 0.003763921 Up hsa_circRNA_103973 hsa_circ_0074269 1.232092034 2.67E-05 0.004725186 Up hsa_circRNA_103384 hsa_circ_0065898 −3.14143122 4.51E-08 0.000135778 Down hsa_circRNA_103677 hsa_circ_0070190 −2.860644018 5.20E-07 0.000391072 Down hsa_circRNA_102031 hsa_circ_0042986 −2.001330988 4.07E-06 0.001200045 Down hsa_circRNA_001259 hsa_circ_0000077 −1.955005875 6.52E-06 0.001634751 Down hsa_circRNA_101308 hsa_circ_0031027 −4.285365461 2.39E-05 0.004486465 Down hsa_circRNA_102050 hsa_circ_0043280 −2.033106203 9.63E-05 0.009577724 Down hsa_circRNA_100319 hsa_circ_0013912 −1.692810166 0.000134465 0.009577724 Down hsa_circRNA_102155 hsa_circ_0045114 −1.06840769 0.000126965 0.009577724 Down hsa_circRNA_100748 hsa_circ_0020926 −1.853767261 0.000223414 0.012679787 Down hsa_circRNA_000593 hsa_circ_0000550 −1.283622936 0.000240927 0.012941213 Down [72]Open in a new tab Abbreviations: DEcircRNAs, differentially expressed circRNAs; FC, fold change; FDR, false discovery rate; CSCC, cervical squamous cell cancer. Functional Annotation of Host Genes of DEcircRNAs A total of 73 host genes of DEcircRNAs were obtained. GO analysis indicated that these host genes were significantly enriched in platelet activation (FDR = 1.32E-02), protein phosphorylation (FDR = 1.35E-02), nucleus (FDR = 6.49E-06) and protein binding (FDR = 4.48E-04) ([73]Figure 2A–C). According to the KEGG pathway enrichment analysis, several pathways were significantly enriched, including Protein processing in endoplasmic reticulum (FDR = 7.74E-04), RNA polymerase (FDR = 3.15E-02) and RNA transport (FDR = 3.42E-02) ([74]Figure 2D). Figure 2. [75]Figure 2 [76]Open in a new tab Significantly enriched GO terms and KEGG pathways of host genes of DEcircRNAs between CSCC tumour tissues and normal cervix tissues. (A) BP, biological process; (B) CC, cellular component; (C) MF, molecular function; (D) KEGG pathways. The x-axis shows counts of host genes enriched in GO terms or KEGG pathways and the y-axis shows GO terms or KEGG pathways. The colour scale represented -lg FDR. CeRNA (DEcircRNA-DEmiRNA-DEmRNA) Regulatory Network The ceRNA network contained 3 circRNA-miRNA pairs and 158 miRNA-mRNA pairs, including 3 circRNAs, 3 miRNAs, and 138 mRNAs ([77]Figure 3). GO enrichment analysis revealed that these DEmRNAs in ceRNA regulatory network were significantly enriched in DNA replication (FDR = 2.84E-14), mitotic cell cycle (FDR = 4.05E-10), nucleus (FDR = 1.50E-25) and protein binding (FDR = 4.52E-14) ([78]Figure 4A–[79]C). According to the KEGG pathway enrichment analysis, several pathways were significantly enriched, including Cell cycle (FDR = 1.23E-07), p53 signaling pathway (FDR = 2.71E-04) and DNA replication (FDR = 4.60E-03) ([80]Figure 4D). Figure 3. [81]Figure 3 [82]Open in a new tab CeRNA (DEcircRNA-DEmiRNA-DEmRNA) regulatory network. The inverted triangles, rectangle and elliptical nodes indicate DEcircRNAs, DEmiRNAs and DEmRNAs, respectively. Red and blue colour represent up-regulation and down-regulation, respectively. Nodes with black border were DEcircRNA/DEmiRNA/DEmRNA derived from top 10 up- and down-regulated DEcircRNA/DEmiRNA/DEmRNA between CSCC tumour tissues and normal cervix tissues. Figure 4. [83]Figure 4 [84]Open in a new tab Significantly enriched GO terms and KEGG pathways of DEmRNAs in ceRNA regulatory network. (A) BP, biological process; (B) CC, cellular component; (C) MF, molecular function; (D) KEGG pathways. The x-axis shows counts of DEmRNAs in ceRNA regulatory network enriched in GO terms or KEGG pathways and the y-axis shows GO terms or KEGG pathways. The colour scale represented -lg FDR. Validation in qRT-PCR and TCGA Database Three DEmRNAs (CCNB2, RRM2 and CDKN2A), three DEmiRNAs (hsa-miR-199b-5p, hsa-miR-125b-5p and hsa-let-7c-5p) and one DEcircRNA (hsa-circ-0000069) were selected to perform validation in qRT-PCR and TCGA. Based on our integrated analysis, CCNB2, RRM2, CDKN2A and hsa-circ-0000069 were up-regulated while hsa-miR-199b-5p, hsa-miR-125b-5p and hsa-let-7c-5p were down-regulated in CSCC. The qRT-PCR results were consistent with those of our integrated analysis ([85]Figure 5). In addition, the validation of CDKN2A, CCNB2, hsa-let-7c-5p and hsa-miR-125b-5p in TCGA database was consistent with our integrated analysis as well ([86]Figure S1). There was no data available for hsa_circ_0020594 and hsa_circ_0000069 in the TCGA datasets. Based on the results of survival analysis, the expression of hsa-let-7c-5p (p < 0.0001) and hsa-miR-125b-5p (p = 0.00045) was significantly correlated with the overall survival time of CSCC patients ([87]Figure S2). Figure 5. [88]Figure 5 [89]Open in a new tab The quantitative real-time polymerase chain reaction (qRT-PCR) results of the DEmRNAs, DEmiRNAs and DEcircRNAs in CSCC tumour tissues. The x-axis represents the DEmRNAs/DEmiRNAs/DEcircRNAs and the y-axis represents log[2] (fold change). The p-value was assessed by one-way analysis of variance (ANOVA). *Indicates p < 0.05 and **Indicates that p < 0.01. Discussion CSCC, with a high recurrence rate and poor prognosis, seriously affects women’s health worldwide.[90]^10 Despite the rapid development of study on CSCC, the five-year survival rate remains not optimistic.[91]^11 As the specific molecular mechanism remains unclear, it is of great clinical diagnostic significance for illuminating the specific role of circRNAs by serving as ceRNAs in CSCC. CDKN2A (full name: cyclin-dependent kinase inhibitor 2A), mapped to chromosome 9, codes for the p16INK4A and p14ARF proteins, and can induce cell cycle arrest.[92]^12 Liu et al linked overexpression of p16INK4A with poor prognosis of CC.[93]^13 Niu et al found that up-regulation of CDKN2A was observed in both high-grade squamous intraepithelial lesions (HSILs) and squamous cell carcinomas (SCCs), which suggested that CDKN2A may contribute to the development of CC.[94]^2 Our analysis, qRT-PCR validation and TCGA validation results consistently showed higher expression levels of CDKNA2 in CSCC, which indicated its critical role in CSCC. The B-type cyclins, including B1 (CCNB1) and B2 (CCNB2), belong to cyclin family and are essential for cell cycle regulation. Increased CCNB2 is associated with various types of cancer, including CC.[95]^14 Compared with normal and CIN1/CIN2, CCNB2 was identified to be highly expressed in tumours and CIN3/CIS by using a microarray technique, which indicated the crucial role of CCNB2 in the early phase of tumorigenesis.[96]^15 A previous study reported that overexpressed CCNB2, belonging to the mitosis pathway, was involved in cell cycle and closely associated with CC.[97]^14 Similarly, a comparative analysis suggested that up-regulated CCNB2 promoted cell cycle transition in CC as well.[98]^16 In this analysis, up-regulated CCNB2 was detected in bioinformatics analysis, qRT-PCR validation, and TCGA validation results, which suggested the importance of CCNB2 in CSCC. It has been reported that hsa-miR-125b-5p showed controversial properties in different tumours. For instance, compared with adjacent tissues, hsa-miR-125b-5p showed up-regulated trend in stage I lung adenocarcinoma tissues.[99]^17 Down-regulation of hsa-miR-125b-5p was observed in bladder cancer, which may be involved in the development of bladder cancer.[100]^18 Decreased hsa-miR-125b-5p was observed in CC cells.[101]^19 In contrast, hsa-let-7c-5p showed consistent trends in a variety of tumours. Down-regulated hsa-let-7c-5p was observed in endometrioid endometrial carcinoma.[102]^20 Significant decrease hsa-let-7c-5p was detected in prostatic cancer.[103]^21 QRT-PCR analysis showed that hsa-let-7c-5p was down-regulated in grade 3 endometrial cancer.[104]^22 Let-7c showed a down-regulated trend in cervical intraepithelial lesions.[105]^23 In addition, both hsa-miR-125b-5p and hsa-let-7c-5p were indicated to be related to the prognosis of colon cancer.[106]^24 Dysregulated hsa-miR-125b-5p and hsa-let-7c-5p were detected by three approaches, including integrated analysis, validation in qRT-PCR and validation with TCGA, in the current study, which indicated their crucial roles in CSCC. In addition, both hsa-miR-125b-5p and hsa-let-7c-5p were determined to be significantly correlated with the overall survival time of CSCC patients. In 2016, Guo et al reported that overexpressed hsa_circ_0000069 was observed in colorectal cancer and correlated to patients’ age and TNM stage.[107]^25 As expected, knockdown of hsa_circ_0000069 could inhibit cell proliferation, migration, and invasion in vitro.[108]^25 However, there is no previous report on hsa_circ_0020594. In this analysis, CDKN2A was targeted by hsa-miR-125b-5p, and CCNB2 was targeted by hsa_circ_0020594, respectively. Then, we found that hsa_circ_0000069 and hsa_circ_0020594 may act as ceRNAs to capture hsa-miR-125b-5p and hsa-let-7c-5p. Therefore, we concluded that the regulation of the hsa_circ_0000069/hsa-miR-125b-5p/CDKN2A and hsa_circ_0020594/hsa-let-7c-5p/CCNB2 axes may be involved in CSCC. This study laid the foundation for further research on ceRNAs of CSCC, and these findings need to be verified in the future. Our study also had a limitation. The sample size used for qRT-PCR validation in this study was small. Although the results were validated by qRT-PCR and TAGA analysis and shown to generally exhibit the same pattern as those in our integrated analysis, studies with larger sample sizes are needed to confirm these findings. Funding Statement This study was supported by Key Special Funding of Health and Family Planning Commission of Changning District, Shanghai: Specialization on Gynecological Minimal Invasive (20151005) and Scientific Research Project of Shanghai Health Committee (201940252). The funding body had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript. Disclosure The authors report no conflicts of interest in this work. References