The profile analysis of circular RNAs in cervical cancer

Abstract Cervical cancer (CC) is the third most common cancer among women and has a high mortality rate at the advanced stage. The mechanisms underlying the development and progression of CC are still elusive. Circular RNAs (circRNAs) play an important role in various physiological and pathological processes. The aim of this study was to identify the circRNAs significantly associated with cervical squamous cell carcinoma (CSCC), in order to discover novel diagnostic markers and elucidate their mechanistic basis. The circRNA expression profiles of CSCC and paired para-cancerous cervical tissues was downloaded from the Gene Expression Omnibus. Bioinformatics analysis were used to screen for the differentially expressed circRNAs (DECRs). The expression levels of hsa_circ_0000745, hsa_circ_0084927, hsa_circ_0002762, hsa_circ_0075341, hsa_circ_0007905, hsa_circ_0031027, hsa_circ_0065898, hsa_circ_0070190, and hsa_circ_0078383 were verified in CC and normal cervical tissues by quantitative real-time PCR. A total of 197 DECRs were identified between the CSCC and normal tissues, including 87 upregulated and 110 downregulated circRNAs. In addition, 37 miRNAs were predicted for the upregulated circRNAs and 39 for the downregulated circRNAs. Functional analysis showed that the DECRs were associated with positive regulation of substrate adhesion-dependent cell spreading, metabolism, positive regulation of GTPase activity, protein regulation, and intercellular adhesion. The MAPK signaling pathway that plays a significant role in the progression of CC, was also enriched. Consistent with the in-silico analysis, hsa_circ_0000745, hsa_circ_0084927, hsa_circ_0002762, hsa_circ_0007905 were upregulated and hsa_circ_0078383 was downregulated in CC tissues (P < .001), whereas hsa_circ_0075341 (P < .001) and hsa_circ_0031027 (P = .001) showed opposite trends. We identified novel diagnostic and therapeutic biomarkers of CSCC along with the mechanistic basis.


Introduction
Cervical cancer (CC) is the third most common cancer among women worldwide, with over 569,800 new cases and 311,400 deaths recorded in 2018 alone. [1] Cervical squamous cell carcinoma (CSCC) accounts for more than 85% of all cases. Although HPV vaccines have shown encouraging results in preventing CC, it remains the fourth most common cause of cancer-related deaths in women. [1,2] While the 5-year survival rate of localized CC is 91.5%, it declines to 16.5% following metastasis. [3,4] The satisfactory prognosis in the early stages of CC can be attributed to the recent advances in conventional treatment strategies. [5][6][7] However, approximately two-thirds of the CC patients are diagnosed at the advanced stage due to lack of practical diagnostic biomarkers, [8] which makes the standard treatment strategies ineffectual. Therefore, it is essential to screen for new diagnostic/therapeutic biomarkers in order to improve prognosis of CSCC.
Circular RNAs (circRNAs) are a category of noncoding RNAs formed via the head-to-tail splicing of exons, which result in covalently closed continuous loops lacking both the 5 0 caps and 3 0 polyadenylated tails. [9,10] Recent studies show that circRNAs regulate tumor cell proliferation, survival and metastasis as tumor suppressors or oncogenes. [11] Hansen et al [12] showed that circRNAs, acting as microRNA sponges and regulate the expression of target mRNAs post-transcriptionally. The circRNA BCRC-3 inhibits bladder cancer proliferation via the miR-182-5p/ p27 axis, [13] while circABCB10 acts as an oncogene in osteosarcoma by sponging miR-203, and facilitates nonsmall cell lung cancer cell progression and migration through the miR-1252/FOXR2 axis. [14,15] Li et al [16] recently reported that several circRNAs are aberrantly expressed in the early stage of bladder cancer, and are therefore promising biomarkers for early diagnosis. In addition, there is also evidence indicating both oncogenic and suppressive functions of specific circRNAs in CC. For instance, circ0023404 inhibits CC cell proliferation, cell cycle progression, migration, and invasion by targeting the mir-136/TFCP2/YAP pathway, [17] while circRNA8924, promotes the malignant and invasive behavior of CC cells by modulating CBX8 as a competitive endogenous (ce) RNA to miR-518d-5p/519-5p. [18] Taken together, it is essential to screen for the circRNAs associated with CC development and progression, in order to identify potential biomarkers for early diagnosis and targeted therapy. Li et al [19] recently profiled circRNAs using the microarray tool. To this end, we obtained the circRNA expression data of CSCC and healthy tissues from the Gene Expression Omnibus (GEO) database, and identified the differentially expressed circRNAs (DECRs) in CSCC. In addition, we also identified the functions of the relevant circRNAs, and elucidated the circRNA-miRNA network potentially involved in CSCC progression.

Identification of DECRs and the target miRNAs
The circRNA expression data of 5 pairs of CSCC and the corresponding para-carcinoma cervical tissues, based on the GPL19978 Agilent-069978 Arraystar Human CircRNA microarray V1 platform (Agilent Technologies Inc., MD), was downloaded from the GEO database (accession no. GSE102686).GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/) of the R limma package was used to screen for the DECRs in CSCC, with log jFCj > 1 and P < .05 as the thresholds. [20,21] The target miRNAs of the DECRs were predicted using the online tool CircInteractome (https://circinteractome.nia.nih.gov/). [22] In addition, studies published till January 20, 2019 on CC-related miRNAs were extracted from Pubmed using the key words "cervical cancer" and "miRNA". Finally, the co-expressing miRNAs of the 2 datasets were used as the target miRNAs of DECRs.

Construction of the circRNA-miRNA-target gene interaction network
The circRNA-miRNA-target gene interaction networks were constructed using the Cytoscape software (3.6.1 version), [23] with the top 5 up-and down-regulated circRNAs and their putative target miRNAs and mRNAs. The miRNAs were predicted using CircInteractome, and the 5 highest scoring (highest context + score percentile) miRNAs were included in the network. In addition, the co-expressing miRNAs extracted from the literature search were also included. The mRNAs targeted by the miRNAs were predicted using miRanda (http://www.microrna.org/) and the top 5 mRNAs (with the lowest mirSVR score) were selected, along with those obtained from literature search. To further validate the analysis, some target mRNAs were verified using the Oncomine platform.

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses
The Database for Annotation, Visualization, and Integrated Discovery (https://david.ncifcrf.gov/, version 6.8) online tool [24,25] was used to identify the origin genes of the DECRs, and their putative functions were determined using the gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. The gene functions were determined on the basis of the 3 GO termsbiological process, cellular component (CC), and molecular function. The associated pathways were elucidated using KEGG analysis. P < .05 indicated the significantly enriched GO terms and KEGG pathways.

Samples collection
A total of 22 CSCC and 9 normal cervical tissue samples were collected from August 2018 to September 2019 at the First Affiliated Hospital of Guangxi Medical University for quantitative real-time PCR (qRT-PCR). The CSCC patients were in stages IB1-IIA1, did not present any other malignancies, and had not received chemotherapy or radiotherapy before samples collection. The diagnosis of CSCC was affirmed by a pathologist, and written informed consent was obtained from all patients. The study was approved by the Ethics Committee of the first affiliated Hospital of Guangxi Medical University.

RNA isolation and reverse transcription
Total RNA was extracted from the tissues using the TRIzolReagen (TAKARA) according to the manufacturer's instructions, and quantified using NanoDropND-2000 Lite spectrophotometer (Thermo Scientific). The purity and concentration of each RNA sample were in the normal range. CircRNA cDNA was synthesized using the PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Dalian, Liaoning, China).

Quantitative real-time PCR (qRT-PCR)
QRT-PCR was performed using the TB Green Premix Ex Taq II (TaKaRa) on the Applied Biosystems7500 Fast Real-Time PCR System (Thermo Fisher Scientific). The reaction conditions were 95°C for 30 seconds, and 40 cycles at 95°C for 5 seconds and 60°C for 30 seconds. Melting curves were plotted and the relative expression levels were calculated using the 2 ÀDDCT method with  Table 1.

Statistical analysis
Statistical analysis was performed using SPSS 20.0 software (IBM, Chicago, IL). Independent t test or Chi-square test was used to compare different groups. P value < .05 was considered statistically significant.

CSCC is associated with an altered circRNA profile
We identified a total of 197 DECRs in the CSCC tissues relative to normal cervix tissues, of which 87 were up-regulated and 110 down-regulated. The DECRs were validated on the basis of their sequences using CircBase (http://www.circbase.org/). The top 40 DECRs are listed in Table 2, and their chromosomal distribution is illustrated in Figure 1. To further determine the trend of these DECRs in CSCC, we performed hierarchical clustering analysis using R package "gplot", and the volcano plot is shown in Table 1 Primer sequences of the circRNAs.

Functional analysis of DECRs
The DECRs are derived from 182 genes, and their functional analysis revealed the following most significantly enriched GO terms: molecular functionprotein binding, GTPase activator activity and cadherin binding involved in cell-cell adhesion, biological processpositive regulation of substrate adhesiondependent cell spreading, cell-cell adhesion and positive regulation of GTPase activity, and CCmembrane, nucleoplasm, and cytosol (Fig. 3). The KEGG pathway analysis indicated that endocytosis, focal adhesion, and protein processing in endoplas-mic reticulum, and the MAPK signaling pathway were the most enriched. The latter is significantly correlated with CSCC, and includes the CDC42, MAPKAPK5, MAPK9, FLNB, NFATC3, FLNA, and PRKCB genes (Fig. 4). We next constructed a circRNA-miRNA-mRNA interaction network to better elucidate the function of the DECRs in CSCC. The ceRNA network of the 5 most significant up-and down-regulated circRNAs is shown in Figure 5.

Validation of DECRs
To verify the results of circRNAs sequencing and determine the possible mechanisms of circRNAs in the development of the occurrence of CSCC, the expression levels of 5 up-regulated and 4 down-regulated DECRs were verified. A total of 22 CSCC and 9 normal cervical tissue samples were used for validation, and each sample was repeated 3 times PCR. Hsa_circ_0000745, hsa_circ_0084927, hsa_circ_0002762, hsa_circ_0007905 were highly upregulated ( Fig. 6A-C, E) and hsa_circ_0078383 was     Fig. 2I) in CSCC tissues, which was consistent with the high-throughput sequencing results. In contrast, the expression levels of hsa_circ_0065898, and hsa_circ_0070190 were similar between the CSCC and normal tissues (P = .48, Fig. 6G; P = .31, Fig. 6H). In addition, hsa_circ_0075341 (P < .05, Fig. 6D) and hsa_circ_0031027, (P < .05, Fig. 6F) were inversely expressed in CSCC and normal tissues.

Discussion
The rapid development of bioinformatics and high-throughput sequencing in recent years has helped identify several circRNAs, which regulate gene expression by acting as ceRNAs that "absorb" target miRNAs. [12] Studies have shown a clear role of miRNAs in CC progression and prognosis. For instance, Zhou et al [26] revealed that miRNA-218 inhibited CC by blocking immune escape of the tumor cellsvia IDO1 downregulation. Persistent infection of the cervical squamous or glandular epithelium by oncogenic human papillomaviruses is a major risk factor of CC. Zheng et al [27] showed that HPV16 E7 regulated the expression profiles of circRNAs in CC cells. Although circRNAs are increasingly being considered as biomarkers for tumor diagnosis and targeted therapy, [28,29] not much is known regarding their role in CC. To this end, we identified 197 DECRs between CSCC and para-cancer tissue with GEO datasets, and further validated the overexpression of hsa_circ_0000745, hsa_circ_0084927, hsa_circ_0002762, and hsa_circ_0007905, and the downregulation of hsa_-circ_0078383 in CSCC tissues relative to normal tissues. In contrast, the expression of hsa_circ_0065898, hsa_-circ_0070190, hsa_circ_0075341, and hsa_circ_0031027 in CSCC and normal cervical tissues need further verification with a larger sample size. At the same time, the errors of highthroughput sequencing cannot be ruled out. Eighty-seven circRNAs were significantly upregulated in CSCC, including hsa_circ_0000745, has_circ_0023404, and hsa_circ_0084927. Hsa_circ_0000745 was targeted by hsa-miR-136, hsa-miR-1229, hsa-miR-1273, and hsa-miR-618. Through Oncomine analysis, we found that the target genes of hsa-miR-136 included E2F1, TFCP2 and RBM41, all of which are overexpressed in CSCC. Lu et al [30] reported that miR-136 inhibited proliferation and promoted apoptosis and radio-sensitivity in cervical carcinoma cells by targeting E2F1 through the NF-kB signaling pathway, and was associated with improved prognosis. Recently, Jiao et al [31] discovered that hsa_circ_0000745 promotes the development of CC. In addition, hsa_circ_0023404 inhibits proliferation, cell cycle progression, migration, and invasion of CC cells via the mir-136/TFCP2/YAP pathway. [17] Hsa_circ_0084927 was targeted by hsa-miR-520h, hsa-miR-874, hsa-miR-634, hsa-miR-1250, and hsa-miR-520g. Studies have correlated mir-520h to the occurrence, progression, invasion, and metastasis of ovarian, breast, lung, and pancreatic cancers. [32][33][34][35] The target genes of hsa-miR-520h were KLHL2, ATPBD4, LHX8, and ZNF318, which are also expressed at high levels in CSCC. Chang et al [36] showed that overexpression of miR-520h is related to CC metastasis. Taken together, overexpression of hsa_circ_0000745, hsa_circ_0023404, and hsa_-circ_0084927 promotes CSCC cell differentiation, proliferation, and metastasis. Hsa_circ_0078383 was targeted by hsa-miR-520h, hsa-miR-1197, hsa-miR-1298, and hsa-miR 520g. The target mRNAs of hsa-miR-520h included PP2A/C, KLHL2, ATPBD4, LOC647979, LHX8, ZNF318, and SMOC2, of which PP2A, KLHL2, ATPBD4, ZNF318, and SMOC2 are significantly downregulated in CC by verified in Oncomine platform. Chang et al. [36] showed that miR-520h promotes metastasis of CC cells by downregulating PP2A. Hsa_circ_0065898 was significantly downregulated in the CSCC tissues, and was targeted by hsa-miR-1245, hsa-miR-1252, hsa-miR-146b-3p, hsa-miR-873, and hsa-miR-874. The target mRNAs of hsa-miR-146b-3p included RUNDC2B, RUNDC2C, PAK7, ANKRD11, LOC203274, and HPGD, all of which except LOC203274, are significantly downregulated in CC. Yao et al [37] reported that miR-146b-3p promotes the development of CC by downregulating HPGD. The parent genes of the DECRs were significantly enriched for functions like positive regulation of substrate adhesiondependent cell spreading, metabolism, positive regulation of GTPase activity, protein regulation, and intercellular adhesion, all of which are related to the proliferation and invasion of tumor cells. KEGG pathway analysis showed significant enrichment of the mitogen-activated protein kinases (MAPK) signaling pathway. Dysregulated MAPK signaling pathway is associated with cancer-related functions like tumor cell proliferation, differentiation, migration, senescence, and apoptosis. [38] PRKCB, the parent gene of hsa_circ_0038645, regulates the MAPK signaling pathway by phosphorylating the RAS signaling pathway (Fig. 4). The expression of PRKCB was related to the survival of CC patients in TGCA database (P = .004). Since both pathways influence and regulate each other, we surmise that the Ras-MAPK axis plays a key role in CSCC progress and metastasis.

Conclusions
We identified 197 DECRs in CSCC, and established the circRNAassociated ceRNA network that can help unearth novel diagnostic and therapeutic biomarkers, and elucidate the mechanistic basis of CSCC.