Reconstruction and analysis of circRNA‑miRNA‑mRNA network in the pathology of cervical cancer.

The present study was performed with the aim of understanding the mechanisms of pathogenesis and providing novel biomarkers for cervical cancer by constructing a regulatory circular (circ)RNA‑micro (mi)RNA‑mRNA network. Using an adjusted P-value of <0.05 and an absolute log value of fold-change >1, 16 and 156 miRNAs from GSE30656 and The Cancer Genome Atlas (TCGA), 5,321 mRNAs from GSE63514, 4,076 mRNAs from cervical squamous cell carcinoma and endocervical adenocarcinoma (from TCGA) and 75 circRNAs from GSE102686 were obtained. Using RNAhybrid, Venn and UpSetR plot, 12 circRNA‑miRNA pairs and 266 miRNA‑mRNA pairs were obtained. Once these pairs were combined, a circRNA‑miRNA‑mRNA network with 11 circRNA nodes, 4 miRNA nodes, 153 mRNA nodes and 203 edges was constructed. By constructing the protein‑protein interaction network using Molecular Complex Detection scores >5 and >5 nodes, 7 hubgenes (RRM2, CEP55, CHEK1, KIF23, RACGAP1, ATAD2 and KIF11) were identified. By mapping the 7 hubgenes into the preliminary circRNA‑miRNA‑mRNA network, a circRNA‑miRNA‑hubgenes network consisting of 5 circRNAs (hsa_circRNA_000596, hsa_circRNA_104315, hsa_circRNA_400068, hsa_circRNA_101958 and hsa_circRNA_103519), 2 mRNAs (hsa‑miR‑15b and hsa‑miR‑106b) and 7 mRNAs (RRM2, CEP55, CHEK1, KIF23, RACGAP1, ATAD2 and KIF11) was constructed. There were 22 circRNA‑miRNA‑mRNA regulatory axes identified in the subnetwork. By analyzing the overall survival for the 7 hubgenes using the Gene Expression Profiling Interactive Analysis tool, higher expression of RRM2 was demonstrated to be associated with a significantly poorer overall survival. PharmGkb analysis identified single nucleotide polymorphisms (SNPs) of rs5030743 and rs1130609 of RRM2, which can be treated with cladribine and cytarabine. RRM2 was also indicated to be involved in the gemcitabine pathway. The 5 circRNAs (hsa_circRNA_000596, hsa_circRNA_104315, hsa_circRNA_400068, hsa_circRNA_101958 and hsa_circRNA_103519) may function as competing endogenous RNAs and serve critical roles in cervical cancer. In addition, cytarabine may produce similar effects to gemcitabine and may be an optional chemotherapeutic drug for treating cervical cancer by targeting rs5030743 and rs1130609 or other similar SNPs. However, the specific mechanism of action should be confirmed by further study.

Introduction

Cervical cancer is the fourth most common cause of fatality that threatens women's health excessively in developing counties (1). Cervical cancer cases were associated with 265,700 fatalities in 2012 (2). Approximately half a million new cases occur worldwide each year (2). Although 80% of early-stage cervical cancer cases can be treated with surgery, radiotherapy or chemotherapy, there is still a large number of advanced-stage patients with poor prognoses (3). The association between human papillomaviruses (HPV) and cervical cancer has been strongly verified (4); however, the specific mechanisms remain to be fully elucidated. Therefore, further studies on the underlying mechanism of tumor initiation and development are necessary.

Non-coding (nc)RNA is a functional RNA molecule that is transcribed from DNA but not translated into proteins (5). Micro(mi)RNAs are a vital component of the endogenous ncRNA family, which are ~19–25 nucleotides in length. The majority of miRNAs are highly conserved in sequence and are involved in multiple cellular functions via the post-transcriptional regulation of gene transcription (6). According to previous studies, several miRNAs (7–11) have been consistently reported to be involved in the development of cervical cancer, which suggests that miRNAs are highly correlated with the pathogenesis of cervical cancer.

Competing endogenous RNAs (ceRNAs) are transcripts with the same miRNA response element (MRE) that act as miRNA sequestering molecules and compete to bind to miRNAs to regulate their target genes, thereby affecting the biological behavior of tumors (12). Circular (circ)RNAs are a novel class of the ceRNA and a distinctive type of ncRNA that were identified in plant viruses in 1976 (13). Unlike traditional linear RNA, circRNAs are structured as a distinct closed loop. Therefore, circRNAs are not easily degraded and stably expressed in various organisms (14). Following the binding to miRNA by MREs, circRNAs can also serve an important role in the occurrence and development of tumors as tumor suppressors or proto-oncogenes (15). Recently, various studies have indicated that circRNAs are involved in the initiation and development of multiple diseases, including gastric cancer (16), osteoarthritis (17) and hepatocellular carcinoma (18).

As the regulatory gene for miRNAs, more studies have suggested that circRNAs may serve key roles in the carcinogenesis of multiple types of cancer (16,18,19). However, research regarding cervical cancer is limited. Wang et al (20) revealed that hsa-circ-0101996 combined with hsa-circ-0101119 in peripheral whole blood was identified as the potential biomarkers for human cervical squamous cell cancer. Gao et al (21) reported that hsa-circ-0018289 was upregulated in cervical cancer and promotes proliferation, migration and invasion of tumor cells. Furthermore, Ma et al (22) indicated that activated has-circ-000284 promotes cell proliferation and invasion in cervical cancer. This evidence supports that circRNA is likely to participate in the development of cervical cancer, and probably by indirectly regulating the expression of target genes through affecting miRNAs.

The molecular mechanism underlying the role of ncRNAs in the carcinogenesis and progression of cervical cancer remains unknown. Therefore, in order to further understand the potential role of ncRNAs in cervical cancer, the differentially expressed RNAs (including circRNA, miRNA and mRNA) were identified by microarray and databases and a regulatory circRNA-miRNA-mRNA network was constructed (Fig. 1). The present study may illuminate the underlying mechanisms of cervical cancer pathogenesis and provide novel biomarkers and targets for cervical cancer.

Figure 1.

Flowchart of the construction and clinical evaluation of circRNA-miRNA-mRNA network. TCGA, The Cancer Genome Atlas; CESC, cervical squamous cell carcinoma and endocervical adenocarcinoma; adj.p, adjusted P-value; FC, fold-change; DE, differentially expressed; GO, Gene Ontology; PPI, protein-protein interaction network; GEPIA, the database of Gene Expression Profiling Interactive Analysis; circRNA, circular RNA; miRNA, microRNA.

Materials and methods

Raw data

Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) is an international public repository for high-throughput microarray and sequence-based data (23). GSE102686 circRNA microarray, GSE30656 miRNA profiles and GSE63514 mRNA datasets for cervical cancer were downloaded from GEO, respectively. The fundamental information for these three profiles is summarized in Table I. In addition, the miRNA and mRNA expression datasets for cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), which contained 309 cervical cancer samples with 3 normal samples, were also downloaded from The Cancer Genome Atlas (TCGA) by cBioPortal (http://www.cbioportal.org/) (24).

Table I.

Basic information of the three microarray datasets from Gene Expression Omnibus.

Profile	RNA type	Platform	Organism	Experiment type	Sample size (T/N)	Region	Year
GSE102686	circRNA	GPL19978	Homo sapiens	Non-coding RNA profiling by array	5/5	China	2017
GSE30656	miRNA	GPL6955	Homo sapiens	Non-coding RNA profiling by array	19/10	Netherlands	2012
GSE63514	mRNA	GPL570	Homo sapiens	Expression profiling by array	28/24	USA	2015

T, tumor; N, normal.

Identification of differentially expressed (DE)-miRNAs, DE-circRNAs and DE-mRNAs

The expression difference between the normal and cervical cancer groups was utilized to determine the DE-circRNAs DE-miRNAs and DE-mRNAs. The adjusted P-value (adj.p) and the absolute log value of fold-change (log|FC|) were calculated using R software with limma package (25). The criteria of adj.p <0.05 and log|FC| >1 were adopted to select the DE-circRNAs, DE-miRNAs and DE-mRNAs.

Construction of miRNA-mRNA pairs

Targeted mRNAs of the collected miRNAs were predicted using miRWalk (Version 2.0; http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/) (26), which included 12 databases (Microt4, miRWalk, mir-bridge, miRanda, miRDB, miRMap, Pictar2, PITA, miRNAMap, RNAhybrid, RNA22, and TargetScan). To increase the accuracy of prediction, the targeted genes were select using the criteria as follows: i) The targeted genes should be predicted by TargetScan (27) and miRanda (28); and ii) the targeted genes should be overlapped in at least 8/12 databases. The selected targeted mRNAs were merged with DE-mRNAs of CESC and GSE63514. The overlapped gene sets were analyzed with UpSetR (29) and Venn Plot, and the pairs of miRNAs and mRNAs were subsequently constructed.

Construction of circRNA-miRNA pairs

To further predict the target circRNAs, the sequences of the candidate DE-miRNAs and DE-circRNAs were downloaded from miRbase (30) and GSE102686, respectively. Following this, the miRNA targets of circRNAs were predicted, and the minimum free energy of circRNA-miRNA duplexes was calculated using the RNAhybrid program (https://bibiserv2.cebitec.uni-bielefeld.de/rnahybrid) (31). miRNA target binding sites on the whole circRNA sequences were predicted. To obtain high-quality circRNAs acting as miRNA targets and distinguish those circRNAs acting as miRNA decoys, the circRNAs that had perfect nucleotide pairing between the 2nd and 8th positions of the 5′ end of miRNA sequences were selected (32) and the circRNA-miRNA pairs were identified simultaneously.

Reconstruction of the circRNA-miRNA-mRNA network

The preliminary circRNA-miRNA-mRNA network was reconstructed by combining the pairs of miRNA-mRNA and miRNA-circRNA. The nodes that could not form a circRNA-miRNA-mRNA axis were removed. The circRNA-miRNA-mRNA network was visualized using Cytoscape software (version 3.6.1) (33).

Gene Ontology (GO) and pathway enrichment analysis

GO and pathway enrichment analyses were performed using The Database for Annotation, Visualization and Integrated Discovery (https://david.ncifcrf.gov/) (34) on the differentially expressed genes (DEGs) in the preliminary circRNA-miRNA-mRNA network. The significantly enriched biological items for biological process (BP), cellular components (CC), and molecular functions (MF) were identified with P<0.05. Pathways with P<0.05 were considered as the significantly enriched pathways. The significant GO items and pathways were visualized using Goplot package (35).

Protein-protein interaction (PPI) network analysis

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (https://string-db.org/) (36) was used to predict the association between target genes in regulatory network analysis. To obtain more accurate results, nodes with interaction score <0.7 and that were not connected to the major network were removed. The Molecular Complex Detection (MCODE) (37) plug-in in Cytoscape was used to analyze PPI network modules, and MCODE scores >5 with >5 nodes were set as cut-off criteria with the default parameters (Degree cut-off ≥2, node score cut-off ≥0.2, K-core ≥2 and max depth = 100). The genes in the cluster were considered hubgenes.

Reconstruction of the circRNA-miRNA-hubgene network

In order to identify the association between DE-circRNAs, DE-miRNAs and hubgenes, the aforementioned hubgenes were mapped into the preliminary circRNA-miRNA-mRNA network and the relevant DE-circRNAs and DE-miRNAs were also extracted. Subsequently, the subnetwork that was identified as the circRNA-miRNA-hubgene network was identified. The structure of circRNAs and the secondary stem-loop structure of miRNAs were also predicted by the cancer-specific circRNA database (CSCD) (38) and Vienna RNA (39), respectively.

Evaluation of overall survival (OS) for hubgenes

To identify the effect of the hubgenes on survival, the Gene Expression Profiling Interactive Analysis (GEPIA) database (40) was utilized to explore the association between these hubgenes and OS. The genes with P<0.05 were considered as critical genes.

Pharmacogenomics analysis for critical genes

Critical genes were analyzed using the database of PharmGkb for the potential single nucleotide polymorphisms (SNPs) and applicable medications, as well as the pharmacogenomics pathway. Specifically, information on clinical and SNP annotations for the critical genes were extracted, including genes names, SNPs, associated medications, efficiency, significance, P-values, genotype association with medications, references and pathways.

Results

Identification of DE-circRNAs, DE-miRNAs and DE-mRNAs

With the criteria of adj.p <0.05 and log|FC| >1 a DE-circRNA dataset consisting of 75 DE-circRNAs, with 43 upregulated and 32 downregulated circRNAs, was identified from GSE102686 (Fig. 2A). A DE-miRNAs dataset, which consisted of 16 and 156 DE-miRNAs, was extracted from GSE30656 and TCGA. By merging these two screening results, 5 miRNAs (hsa-miR-21, hsa-miR-99a, hsa-miR-106b, hsa-miR-15b and hsa-miR-203) were obtained. hsa-miR-106b, hsa-miR-15b and hsa-miR-21 were upregulated whereas hsa-miR-99a was downregulated in CESC and GSE30656. However, hsa-miR-203 was upregulated in CESC but downregulated in GSE30656 (Fig. 2B and C). A DE-mRNA dataset, which consisted of 5,321 and 4,076 DE-mRNAs from GSE63514 and CESC, was determined (Fig. 2D and E). In addition, 3,400 targeted genes predicted by miRWalk with 12 databases were identified. The expression of the 75 DE-circRNAs is indicated in Fig. 3A, and the expression of 5 miRNAs in GSE30656 is indicated in Fig. 3B. After integrating the DE-mRNAs from GSE63514 and CESC with the targeted genes from miRWalk, 155 DE-mRNAs were identified (Fig. 4A). The expression and oncoprint of these 155 genes in GSE63514 and TCGA is indicated in Fig. 4B.

Figure 2.

Volcano plots for DE-circRNAs, DE-miRNAs and DE-mRNAs identified from GEO and TCGA. (A) DE-circRNAs from GSE102686; (B) DE-miRNAs from GSE306565; (C) DE-miRNAs from TCGA. (D) DE-mRNAs from GSE63514; and (E) DE-mRNA from TCGA. DE, differentially expressed; circRNA, circular RNA; miRNA, microRNA; TCGA, The Cancer Genome Atlas; GEO, Gene Expression Omnibus.

Figure 3.

Expression of 75 DE-circRNAs and 5 DE-miRNAs in GEO. (A) DE-circRNAs from GSE102686. (B) DE-miRNAs from GSE306565. DE, differentially expressed; circRNA, circular RNA; miRNA, microRNA; GEO, Gene Expression Omnibus.

Figure 4.

The 155 overlapped genes identified from GSE63514, TCGA and miRWalk with their expression and oncoprint. (A) UpsetR for miRWalk and Venn Plot for the 155 overlapped genes; and (B) expression and oncoprint of the 155 overlapped genes in GSE63514 and TCGA. TCGA, The Cancer Genome Atlas.

By merging the DE-mRNAs of GSE63514 and TCGA with the targeted mRNAs of 5 DE-miRNAs from miRWalk, 266 miRNA-mRNA pairs were selected. Specifically, 72 mRNAs of hsa-miR-15b, 16 mRNAs of hsa-miR-21, 2 mRNAs of hsa-miR-99a, 78 mRNAs of hsa-miR-106b and 58 mRNAs of hsa-miR-203 were obtained.

Using RNAhybrid with the criterion that the circRNAs had perfect nucleotide pairing between the 2nd and 8th positions of the 5′ end of miRNA sequences, 12 circRNA-miRNA pairs were obtained.

As indicated in Fig. 5A, by combining the pairs of miRNA-mRNA and miRNA-circRNA, a preliminary circRNA-miRNA-mRNA network was constructed. The network was composed of 11 circRNA nodes, 4 miRNA nodes, 153 mRNA nodes and 203 edges. The network presented an initial perception of the association between the 11 DECs (hsa_circRNA_103102, hsa_circRNA_002172, hsa_circRNA_101835, hsa_circRNA_102155, hsa_circRNA_102050, hsa_circRNA_103384, hsa_circRNA_104315, hsa_circRNA_400068, hsa_circRNA_000596, hsa_circRNA_101958 and hsa_circRNA_103519), the 4 miRNAs (hsa-miR-21, hsa-miR-203, hsa-miR-15b and hsa-miR-106b) and the 153 mRNAs.

Figure 5.

Preliminary circRNA-miRNA-mRNA, protein-protein interaction and circRNA-miRNA-hubgenes network for the 11 DECs (hsa_circRNA_103102, hsa_circRNA_002172, hsa_circRNA_101835, hsa_circRNA_102155, hsa_circRNA_102050, hsa_circRNA_103384, hsa_circRNA_104315, hsa_circRNA_400068, hsa_circRNA_000596, hsa_circRNA_101958 and hsa_circRNA_103519), the 4 miRNAs (hsa-miR-21, hsa-miR-203, hsa-miR-15b and hsa-miR-106b) and the 155 mRNAs. (A) Preliminary circRNA-miRNA-mRNA network the 11 DECs, the 4 miRNAs and the 155 mRNAs. (B) Protein-protein interaction network for the 155 overlapped genes. Each ellipses represent the cluster identified by MCODE algorithm. (C) The circRNA-miRNA-hubgenes network for the 5 circRNAs (hsa_circRNA_000596, hsa_circRNA_104315, hsa_circRNA_400068, hsa_circRNA_101958 and hsa_circRNA_103519), the 2 mRNAs (hsa-miR-15b and hsa-miR-106b) and the 7 mRNAs (RRM2, CEP55, CHEK1, KIF23, RACGAP1, ATAD2 and KIF11). circRNA, circular RNA; miRNA/miR, microRNA; MCODE, Molecular Complex Detection.

PPI network analysis

The STRING database was used to construct a PPI network based on the 155 overlapped mRNAs. The original network contained 47 nodes and 72 edges. By utilizing the algorithm of MCODE, four clusters were identified (Fig. 5B). Using the criteria of MCODE scores >5 and >5 nodes, one cluster was selected. There were 7 genes (RRM2, CEP55, CHEK1, KIF23, RACGAP1, ATAD2 and KIF11) in this cluster, and these 7 genes were identified as hubgenes.

Construction of the circRNA-miRNA-hubgenes network

By mapping the 7 hubgenes into the preliminary circRNA-miRNA-mRNA network, and also extracting relevant circRNAs and miRNAs, a subnetwork considered as circRNA-miRNA-hubgenes was constructed. This subnetwork consisted of 5 circRNAs (hsa_circRNA_000596, hsa_circRNA_104315, hsa_circRNA_400068, hsa_circRNA_101958 and hsa_circRNA_103519), 2 mRNAs (hsa-miR-15b and hsa-miR-106b) and 7 mRNAs (RRM2, CEP55, CHEK1, KIF23, RACGAP1, ATAD2 and KIF11) with 14 edges (Fig. 5C). There were 22 circRNA-miRNA-mRNA regulatory axes identified in the subnetwork (Table II). The structure and the primary characteristics of the circRNAs in the subnetwork were indicated in Fig. 6 and Table III. In addition, the stem-loop structure of hsa-miR-15b and hsa-miR-106b was revealed (Fig. 7A). The binding site to the circRNAs was also indicated (Fig. 7B).

Table II.

Regulatory axis identified from the circRNA-miRNA-hubgene network.

Regulatory axis	circRNA	miRNA	Hubgene
1	hsa_circRNA_000596	hsa-miR-15b	CHEK1
2	hsa_circRNA_000596	hsa-miR-15b	CEP55
3	hsa_circRNA_104315	hsa-miR-15b	CHEK1
4	hsa_circRNA_104315	hsa-miR-15b	CEP55
5	hsa_circRNA_400068	hsa-miR-15b	CHEK1
6	hsa_circRNA_400068	hsa-miR-15b	CEP55
7	hsa_circRNA_000596	hsa-miR-15b	KIF23
8	hsa_circRNA_000596	hsa-miR-15b	RACGAP1
9	hsa_circRNA_104315	hsa-miR-15b	KIF23
10	hsa_circRNA_104315	hsa-miR-15b	RACGAP1
11	hsa_circRNA_400068	hsa-miR-15b	KIF23
12	hsa_circRNA_400068	hsa-miR-15b	RACGAP1
13	hsa_circRNA_101958	hsa-miR-106b	KIF23
14	hsa_circRNA_101958	hsa-miR-106b	RACGAP1
15	hsa_circRNA_103519	hsa-miR-106b	KIF23
16	hsa_circRNA_103519	hsa-miR-106b	RACGAP1
17	hsa_circRNA_101958	hsa-miR-106b	RRM2
18	hsa_circRNA_103519	hsa-miR-106b	RRM2
19	hsa_circRNA_101958	hsa-miR-106b	ATAD2
20	hsa_circRNA_103519	hsa-miR-106b	ATAD2
21	hsa_circRNA_101958	hsa-miR-106b	KIF11
22	hsa_circRNA_103519	hsa-miR-106b	KIF11

miRNA, microRNA; circRNA, circular RNA.

Figure 6.

Structure of the 5 circRNAs predicted by the cancer-specific CircRNA. (A) hsa_circRNA_000596; (B) hsa_circRNA_400068; (C) hsa_circRNA_104315; (D) hsa_circRNA_101958; (E) hsa_circRNA_103519. The different colors in the outer and inner ring represent the different exons and the number/position of MRE, RBP and ORF. MRE, microRNA response element; RBP, RNA binding protein; ORF, open reading frame; circRNA, circular RNA.

Table III.

Primary characteristics of the 5 circRNAs identified from the circRNA-miRNA-hubgene network.

circRNA	Chromosome	Start position	End position	Strand	Location	Gene symbol	Regulation
hsa_circRNA_000596	Chr15	100589061	100739628	+	Exon	ADAMTS17	Upregulated
hsa_circRNA_400068	Chr22	20109884	20114047	+	Exon	RANBP1	Downregulated
hsa_circRNA_104315	Chr7	16298014	16330476	−	Exon	ISPD	Downregulated
hsa_circRNA_101958	Chr17	4192524	4207021	−	Exon	UBE2G1	Downregulated
hsa_circRNA_103519	Chr3	179042885	179051613	+	Exon	ZNF639	Upregulated

miRNA, microRNA; circRNA, circular RNA.

Figure 7.

Predicted stem-loop structure of hsa-miR-15b and has-miR-106b with the binding site to the corresponding circRNAs. (A) Predicted secondary stem-loop structure of hsa-miR-15b and has-miR-106b. The colors in each base pair represent its probabilities. (B) Binding site between hsa-miR-15b, has-miR-106b, hsa_circRNA_000596, hsa_circRNA_104315, hsa_circRNA_400068, hsa_circRNA_101958 and hsa_circRNA_103519. miRNA/miR, microRNA; circRNA, circular RNA.

GO and pathway enrichment analysis

The 155 overlapped genes were utilized for GO and pathway enrichment analysis using DAVID. For GO analysis, when considering BPs, the DEGs were enriched in intracellular signaling cascade, cell cycle and protein amino acid phosphorylation. With regards to CC, the top three enriched items were cell death, neuron projection and cytoskeleton. In terms of MF, transcription factor activity, transcription regulator activity and protein kinase activity were enriched for the first three places. Pathway enrichment analysis indicated that the top three significantly enriched pathways were pathways in cancer, cell cycle and prostate cancer (Fig. 8).

Figure 8.

Top 10 items of biological process, cellular component and molecular function in Gene Ontology and pathway enrichment analysis for the 155 overlapped genes.

Evaluation of OS for hubgenes

GEPIAT was used to assess the OS for the 7 hubgenes. Notably, higher expression of RRM2 revealed a significantly poorer OS (hazard ratio = 1.7, P=0.029). However, no significant effect was indicated for the remaining 6 hubgenes regarding OS (Fig. 9).

Figure 9.

Overall survival for the 7 hubgenes (RRM2, CEP55, CHEK1, KIF23, RACGAP1, ATAD2 and KIF11) from TCGA. TPM, transcripts per million; HR, hazard ratio; TCGA, The Cancer Genome Atlas.

Pharmacogenomics analysis for RRM2

As RRM2 significantly impacted OS, this gene was selected for the analysis of pharmacogenomics. Using PharmGkb, the SNPs of rs5030743 and rs1130609 were identified, which can be treated with cladribine and cytarabine (Table IV). Notably, cladribine is used for hairy cell leukemia, acute myeloid leukemia and chronic lymphocytic leukemia (41–43) while cytarabine is a chemotherapeutic drug utilized for various types of cancer, such as breast and lung cancer (44,45). Cytarabine, similar to gemcitabine, is activated and metabolized by deoxycytidine kinase and cytidine deaminase (46). In addition, their primary metabolite can also incorporate DNA into cells (47). Hence, the gemcitabine pathway was selected to investigate the mechanism of action of cytarabine through RRM2. As indicated in Fig. 10, gemcitabine produces its therapeutic effect by acting on the enzymes of RRM1, RRM2 and RRM2B, while the RRM1, RRM2 and RRM2B genes also strictly regulate the activity of these enzymes.

Table IV.

Information of SNPs and clinical data for RRM2.

Genes	SNP	Drug	Efficiency	Significant	P-value	Association	Reference	Pathways
RRM2	rs5030743	Cladribine; Cytarabine	Effective	Yes	0.013	Genotypes CG + GG is associated with decreased response to cladribine and cytarabine in children with leukemia, myeloid, acute as compared to genotype CC	PMID: 24024897	Thiopurine pathway; Fluoropyrimidine pathway; Gemcitabine pathway
RRM2	rs1130609	Cladribine; Cytarabine	Efficacy	Yes	0.020	Genotypes GT + TT is associated with decreased response to cladribine and cytarabine in children with leukemia, myeloid, acute as compared to genotype GG	PMID: 24024897	Thiopurine pathway; Fluoropyrimidine pathway; Gemcitabine pathway

SNP, single nucleotide polymorphism.

Figure 10.

Gemcitabine pathway extracted from PharmGkb database. Gemcitabine could exert its therapeutic effect by acting on RRM2, RRM1 and RRM2B. dFdu, 2′,2′-difluorodeoxyuridine; dFdUMP, dFdU monophosphate; NDP, nucleoside diphosphates; dNDP, deoxy-ribonucleoside diphosphate; dNTP, deoxy-ribonucleoside triphosphate.

Discussion

During the past few decades, cervical cancer research has developed rapidly, particularly with regard to the discovery of the etiological factors. Although the key etiological role of HPV in the development of cervical cancer has been reported (48), the specific molecular mechanism remains unclear. Recently, more data has suggested that circRNAs participate in various biological processes (49–52). Dysregulated circRNA expression has been observed in the progression of complex diseases, including cervical cancer (16,18,19,21). circRNAs are a novel class of competing endogenous ncRNAs that are widely expressed in eukaryotic cytoplasm (53). Unlike the RNAs formed by linear splicing, circRNAs have a closed loop structure with no 5′ end cap and 3′ end tail. The majority of known circRNAs are backspliced from exons, which are not easily degraded by exonuclease and have the feature of high abundance, structural stability and tissue specificity (54). With the rapid development of high-throughput RNA sequencing technology accompanied with the extensive data analysis by bioinformatics, circRNAs have been identified to possess the function of miRNA sponges (54–56), gene transcription regulation (54,57) and RNA-binding proteins (58–60) as well as their translation (52,61,62). In addition, studies have confirmed that circRNA serves an important role in tumor cells, so may potentially serve as a novel biomarker and therapeutic target for cancer therapy (63–65). However, more circRNAs need to be uncovered. In the present study, the circRNA microarray profile of GSE102686 for cervical cancer was screened to identify DE-circRNAs. Using the criteria of FDR <0.05 and log|FC| >1 calculated by R with limma package, 75 DE-circRNAs were selected for further analysis.

Previous studies have indicated that circRNA is a type of high-efficiency ceRNA (66,67). It can inhibit the binding of miRNAs to target genes and regulate the expression level of target genes by exerting a miRNA sequestering effect (68). To determine whether the above 75 circRNAs function as ceRNAs in cervical cancer, the targeted miRNAs with the sequence of DE-circRNAs and DE-miRNAs were predicted with RNAhybrid. Furthermore, 12 circRNA-miRNA pairs were identified. hsa-miR-106b, hsa-miR-15b and hsa-miR-21 were identified to be upregulated in CESC and GSE30656, and these results were consistent with previous findings (69–73). Notably, hsa-miR-203 was upregulated in CESC but downregulated in GSE30656. Previous studies have suggested that the expression of miR-203 is typically downregulated in cervical cancer tumors and cell lines (74–76). This finding is in line with the expression profile of GSE30656 but contradicts to TCGA. Reshmi and Pillai (77) also identified that miR-203 can be overexpressed or underexpressed in cervical cancer cell lines. Furthermore, Zhao et al (78) indicated that the expression level of miR-203 in serum of patients with cervical cancer was significantly upregulated; however, miR-203 downregulation was correlated with lymph nodes metastasis. Hence, it was possible that upregulation and downregulation of miR-203 could promote the progression of cervical cancer but serves different functions. Upregulation of miR-203 could also trigger the occurrence of cervical cancer whereas downregulation may enhance the metastatic capacity of lymph nodes. However, the specific mechanism requires further investigation.

Once the 155 overlapping genes with the DE-circRNAs and DE-miRNAs were collected, a circRNA-miRNA-mRNA regulatory network was reconstructed. The 11 circRNAs could bind to hsa-miR-21, hsa-miR-203, hsa-miR-106b or hsa-miR-15b as ceRNAs to regulate the expression of 155 genes. The present results provide evidence for the ceRNA regulatory mechanisms of 11 circRNAs in cervical cancer. The GO and pathway enrichment analysis indicated that 155 genes were involved in various important biological functions and metabolic pathways associated with tumors, including ‘Pathways in cancer’, ‘Cell cycle’, ‘Prostate cancer’, ‘Small cell lung cancer’, ‘p53 signaling pathway’ and ‘Pancreatic cancer’. It also indicated that several type of cancer, such as prostate cancer, lung cancer and pancreatic cancer, share the same pathways with cervical cancer, and provided evidence for the research of pan-cancer. To further understand the functional mechanism of the ceRNA-miRNA-mRNA network, the PPI network was constructed and 7 hubgenes (RRM2, CEP55, CHEK1, KIF23, RACGAP1, ATAD2 and KIF11) were identified from the PPI network. The crucial mechanisms of RRM2, CEP55, CHEK1, ATAD2 and KIF11 in cervical cancer have been previously studied (79–82). However, to the best of our knowledge, KIF23 and RACGAP1 have not yet been investigated. The 22 circRNA-miRNA-mRNA axes, which were identified from the circRNA-miRNA-hubgenes, revealed the competing regulatory associations between the 5 circRNAs and the 2 mRNAs with the 7 hubgenes in cervical cancer. However, as the present study was an in silico research, a further experiment of these 22 regulatory axes is required for validation.

The OS for the 7 hubgenes was evaluated in the present study and it was revealed that RRM2 had a significant effect on OS. Subsequently, RRM2 was investigated using the database of PharmGkb. It was revealed that cytarabine, which is similar to gemcitabine, may produce its therapeutic effect by targeting to the SNPs of rs5030743 and rs1130609 of RRM2 through the gemcitabine pathway. To date, gemcitabine has been indicated to be feasible and effective on cervical cancer (83–85). However, cytarabine has not yet been investigated in cervical cancer, nor the rs5030743 and rs1130609 of RRM2. Although cytarabine differs from gemcitabine in several important respects (86), for instance, gemcitabine differs structurally from cytarabine by a fluorine group substituted at position 2′ on the furanose ring (87), gemcitabine (88–90) and cytarabine (91,92) can be a treatment for breast cancer. Hence, cytarabine may produce a similar effect to gemcitabine and could be an optional chemotherapeutic drug for treating cervical cancer by targeting rs5030743 and rs1130609 or other similar SNPs. However, further study is necessary to validate the hypothesis.

To conclude, the present study constructed and analyzed a circRNA-miRNA-mRNA network based on the ceRNA theory via comprehensive bioinformatics analysis, which may provide some evidence to future studies focused on the molecular mechanisms of cervical cancer. The 5 circRNAs (hsa_circRNA_000596, hsa_circRNA_104315, hsa_circRNA_400068, hsa_circRNA_101958 and hsa_circRNA_103519) may function as ceRNAs to serve critical roles in cervical cancer. In addition, cytarabine may produce a similar effect to gemcitabine and could be an optional chemotherapeutic drug for cervical cancer by targeting rs5030743 and rs1130609 or other similar SNPs. However, the specific mechanism of action should be confirmed by further study.