Comprehensive Analysis of circRNA-miRNA-mRNA Network in Cervical Squamous Cell Carcinoma by Integrated Analysis

Introduction

Cervical cancer (CC) is one of the most widespread gynaecological cancers that threatens women’s health in developing countries.¹ 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.² Human papillomaviruses (HPVs) have been widely considered the main cause of CC.³ Cervical squamous cell carcinoma (CSCC) is the most frequent type, accounting for the majority of CCs.⁴ 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.⁵ 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.⁶ The advent of microarray-based technology has facilitated studies on the molecular mechanisms of different tumours.⁷ 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.⁸ Yi et al identified five circRNAs that may play critical roles in CC.¹ Ma et al revealed that circRNA-000284 is involved in cell proliferation and invasion of CC by sponging miR-506.⁹ 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 (Table 1).

Table 1 List of mRNA/miRNA/circRNA Study Samples from GEO

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 (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 (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 (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 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 (Figure 1A and 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 (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 (Figure 1C). Among them, hsa_circRNA_104052 and hsa_circRNA_103384 were the most up-regulated and down-regulated circRNAs, respectively (Table 3).

Table 2 Top 10 Up- and Down-Regulated DEmRNAs and DEmiRNAs in CSCC

Table 3 Top 10 Up- and Down-Regulated DEcircRNAs in CSCC

Figure 1 The heatmap of top 100 up- and down-regulated DEmRNAs (A), DEmiRNAs (B) and DEcircRNAs (C) between CSCC and normal cervix tissues.

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) (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) (Figure 2D).

Figure 2 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 (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) (Figure 4A–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) (Figure 4D).

Figure 3 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 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 (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 (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 (Figure S2).

Figure 5 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 log2 (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.¹⁰ Despite the rapid development of study on CSCC, the five-year survival rate remains not optimistic.¹¹ 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.¹² Liu et al linked overexpression of p16INK4A with poor prognosis of CC.¹³ 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.² 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.¹⁴ 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.¹⁵ A previous study reported that overexpressed CCNB2, belonging to the mitosis pathway, was involved in cell cycle and closely associated with CC.¹⁴ Similarly, a comparative analysis suggested that up-regulated CCNB2 promoted cell cycle transition in CC as well.¹⁶ 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.¹⁷ Down-regulation of hsa-miR-125b-5p was observed in bladder cancer, which may be involved in the development of bladder cancer.¹⁸ Decreased hsa-miR-125b-5p was observed in CC cells.¹⁹ 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.²⁰ Significant decrease hsa-let-7c-5p was detected in prostatic cancer.²¹ QRT-PCR analysis showed that hsa-let-7c-5p was down-regulated in grade 3 endometrial cancer.²² Let-7c showed a down-regulated trend in cervical intraepithelial lesions.²³ In addition, both hsa-miR-125b-5p and hsa-let-7c-5p were indicated to be related to the prognosis of colon cancer.²⁴ 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.²⁵ As expected, knockdown of hsa_circ_0000069 could inhibit cell proliferation, migration, and invasion in vitro.²⁵ 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.