Regulatory Network Analysis Reveals Prognosis-Related CircRNAs in Pancreatic Cancer Cells


 Objective: To construct a key prognosis-related regulatory network and to identify the epigenomic alterations of RNAs that have crucial functions in cancer pathogenesis.Methods: RNA expression profiles of miRNAs, mRNAs and circRNAs were extracted from the NCBI GEO and EBI ArrayExpress databases. The identification of differentially expressed genes was performed by R language. Databases such as starBase and TCGA were used for annotation, visualization, and integrated discovery. GO and KEGG pathway enrichment analyses were performed by DAVID 6.8. The key prognosis-related regulatory network was constructed with Cytoscape software. The involved circRNAs were verified by quantitative real-time PCR (qRT-PCR) in five pancreatic cancer cell lines. Proliferation ability was assessed by the CCK-8 assay, apoptosis was detected by flow cytometry, and migratory and invasive abilities were assessed by Transwell assays.Results: In total, 798 differentially expressed genes, consisting of 85 circRNAs, 19 miRNAs and 694 mRNAs, were obtained. A miRNA interaction network was predicted and used to construct a prognosis-related regulatory network. GO annotation and KEGG pathway analysis revealed important biological processes and pathways in pancreatic cancers. The key regulatory factors in the network were miR-146b-5p and miR-152. Eight circRNAs included in the network were verified, and hsa_circ_0006502 was downregulated in the five tested pancreatic cancer cell lines and functioned as a tumor suppressor, as predicted.Conclusion: In conclusion, a key prognosis-related regulatory network in pancreatic cancers was constructed through bioinformatics analysis of public RNA databases, and miR-146b-5p and miR-152 were the key regulatory factors. Hsa_circ_0006502 was identified as a tumor suppressor that interacted closely with miR-146b-5p and might serve as a potential therapeutic agent.


Introduction
Pancreatic cancer is one of the most fatal malignancies due to the limited availability of effective therapeutic strategies, and its incidence and mortality rates are similar worldwide. Signi cant progress has been made in the overall ve-year survival rates for many cancers over the past 3 decades but not for pancreatic cancer, mostly because many of the patients are not diagnosed early due to the lack of detectable biomarkers and recognizable symptoms. Therefore, it is urgent to reveal what the mechanisms of tumorigenesis and pathological development processes, such as invasion and metastasis, identify key regulators, and simultaneously screen effective diagnostic biomarkers and potential treatment targets for pancreatic cancer.
Abnormal expression pro les of messenger RNA (mRNA), as well as noncoding transcripts in the human transcriptome, are closely related to human diseases, including pancreatic cancer. Currently, more than 75% of the human genome can be transcribed, but only 2% of transcripts encode proteins [1,2]. The function and underlying mechanism of a large number of noncoding RNAs (ncRNAs) remain unclear.
Circular RNA (circRNA) is an endogenous ncRNA molecule characterized by tissue-speci c expression in a variety of biological cells, structural stability, and gene expression regulatory function [3]. The discovery of a large number of circRNAs and the elucidation of their structure and function offered not only a better understanding of the mechanism of pancreatic cancer but also new strategies for its diagnosis and treatment.
Here, the RNA expression pro les of miRNAs, mRNAs and circRNAs were extracted from the NCBI GEO and EBI ArrayExpress databases. The identi cation of differentially expressed genes was performed by R language. Databases such as starBase and The Cancer Genome Atlas (TCGA) were used for annotation, visualization, and integrated discovery. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed by DAVID 6.8. The key prognosis-related regulatory network was constructed with Cytoscape software. The involved circRNAs were veri ed by quantitative real-time PCR (qRT-PCR) in ve pancreatic cancer cell lines and one normal pancreas cell line. In addition, proliferation ability was assessed by the CCK-8 assay, apoptosis was detected by ow cytometry, and migratory and invasive abilities were assessed by Transwell assays. A prognosis-related competing endogenous RNA (ceRNA) network was thus constructed, and the key regulatory circRNAs involved in the promotion of pancreatic cancer were identi ed and veri ed.

Materials And Methods
Preprocessing of expression pro le data RNA expression pro le datasets that contained no less than 10 pairs of human pancreatic cancer and normal control samples were retrieved from NCBI GEO (http://www.ncbi.nlm.nih.gov/geo/) and EBI ArrayExpress (Https://www.ebi.ac.uk/arrayexpress/). A total of 2 circRNA, 5 miRNA and 7 mRNA expression pro le chip datasets were thus obtained, as shown in Table 1.

Screening Of Differentially Expressed Genes
Flow chat of the bioinformatics analysis was shown in Fig. 1. Signi cantly differentially expressed circRNAs in the two circRNA datasets were screened with the Limma package [3], with a p value less than 0.05 and | logFC |> 0.585 (multiplied by 1.5) as the thresholds. All the circRNA names were converted into those commonly used in the circRNA database (Http://www.circbase.org/) [5](see Table 1). The differentially expressed miRNAs and mRNAs were screened by MetaDE (version 1.0.5, https://CRAN.R_project.org/package=MetaDE) [6] in R3.4.1, with tau 2 = 0 and Q Pval > 0.05 as the homogeneity test parameters and FDR < 0.05 as the threshold for signi cant expression differences between genomes. The MetaQC package (version 0.1.13, https://cran.r_project.org/web/pathe ckages/MetaQC/index.html) [7] in R3.4.1 combined with standardized mean rank (SMR) and principal component analysis (PCA) was used as an objective data quality control measure and con rmed that the data in the 5 miRNA datasets and 7 mRNA datasets were well distributed and suitable for further analysis (supplementary Fig. 1). Bidirectional hierarchical clustering of the above differentially expressed RNAs was performed based on Euclidean distance with the Pheatmap package (version 1.0.8, https://cran.r_project.org/package=pheatmap) [8] in R3.4.1.
The data of pancreatic cancer samples with survival information were obtained from the TCGA database (https://gdc_portal.nci.nih.gov/). The miRNAs and mRNAs in the above regulatory relationships that were signi cantly related to prognosis were screened by single-factor Cox regression analysis with the survival package (version: 2.40-1, https://cran.r-project) in the R3.4.1 language .org / package = survival) [19] to construct a ceRNA network related to prognosis. The miRNA-associated pathways involved in the ceRNA network were downloaded from miRWalk 2.0 (http://zmf.umm.uni_heidelberg.de/apps/zmf/mirwalk2/index.html) [20].
Veri cation of the 8 most differentially expressed circRNAs in pancreatic cancer cell lines The pancreatic cancer cell lines AsPC-1, PANC-1, BxPC-3, CFPAC-1, and HPAF-II and the normal human pancreatic cell line HPDE6-C7 were purchased from the Shanghai Cell Bank of the Chinese Academy of Sciences, thawed and cultured in DMEM (Gibco, C11995500BT) or RPMI-1640 medium (Gibco, C11875500BT) containing 10% fetal bovine serum, 1% penicillin/streptomycin (complete medium, Gibco) with or without 1% glutamine and 1% sodium pyruvate (Gibco) according to the manufacturer's protocols. The cells were cultured in the recommended incubator at 37 °C supplied with 5% CO2.

Statistical analysis
All results are presented as the mean ± standard deviation (SD). Statistical analyses were performed with GraphPad Prism 7 (GraphPad Software, San Diego, CA). Student's t test and one-way ANOVA were utilized to analyze signi cant differences, and P < 0.05 and P < 0.01 were considered to represent statistically signi cant and extremely signi cant differences, respectively.

Results
In total, 798 differentially expressed genes were obtained A total of 387 and 616 differentially expressed circRNAs (supplementary table 2) were obtained from the two circRNA datasets GSE69362 and GSE79634, respectively. Figure 2. A shows the corresponding bidirectional hierarchical clustering heat map. A total of 109 circRNAs appeared in both datasets, among which 85 circRNAs had consistent expression changes, 71 were upregulated and 14 were downregulated ( Fig. 2. B, supplementary table 3). Nineteen differentially expressed miRNAs were obtained (supplementary table 4), and they could successfully distinguish tumors from normal samples (Fig. 2. C).
A total of 694 signi cantly differentially expressed mRNAs were identi ed (supplementary table 5). The bidirectional and hierarchical clustering heat map indicated that these genes were expressed in different types of samples in the 7 datasets (the red and green colors indicate differentially expressed genes), and the differences were obviously consistent ( Fig. 2. D, the distribution of red and green spots in each dataset is very consistent). In total, 798 differentially expressed genes were obtained.

Prediction Of A Mirna Interaction Network
The sequences of the 85 differentially expressed circRNAs were downloaded from circBase, while the sequences of the 19 miRNAs were downloaded from miRbase. There were 48 paired circRNA-miRNA interactions predicted by the miRanda program, including 36 circRNAs (7 downregulated and 29 upregulated) and 7 miRNAs (2 downregulated and 5 upregulated). Overall, there were 43 nodes and 48 edges (Fig. 3. A).
The target genes of the above 7 miRNAs were predicted with TargetScan Release 7.1, MicroCosm version 5, miRecords and miRNAMap ( Fig. 3. B), and 196 miRNA-mRNA interaction pairs were identi ed. The target mRNA interactions were also explored with BioGRID version 3.4.153, HPRD Release 9 and STRING version 10.5, and 86 gene interaction pairs appeared in at least two datasets that were included in the miRNA regulatory network (Fig. 3. C). A miRNA-mRNA interaction network was thus constructed, containing a total of 171 nodes with 7 miRNAs (2 downregulated and 5 upregulated) and 164 mRNAs (126 downregulated and 38 upregulated) and 282 connecting edges with 196 miRNA-targeted mRNA connections and 86 mRNA interaction connections (Fig. 3. D).
GO annotations and KEGG pathway analysis of the target genes identi ed 22 signi cantly related GO functional annotations with 9 biological processes, 4 cellular components and 9 molecular functions ( Fig. 3. E) and 11 KEGG pathways, respectively ( Fig. 3. F, supplementary table 6). The target genes were signi cantly involved in biological processes such as apoptosis and programmed cell death and signaling pathways such as tight junctions and the MAPK signaling pathway.

Construction Of A Prognosis-related Regulatory Network
In total, 3 miRNAs and 75 mRNAs were con rmed to be signi cantly related to pancreatic cancer prognosis. MiR-146b-5p and miR-152 were overexpressed, while miR-653 was downregulated in cancer samples. Lower expression of miR-146b-5p and miR-152 predicted a longer survival time, and higher expression of miR-653 was associated with a good prognosis (Fig. 4. A). A ceRNA network related to prognosis was constructed with the above differentially expressed circRNAs, miRNAs and mRNAs (Fig. 4. B).
Overall, 161 KEGG pathways were obtained from miRWalk 2.0 in which miR-146b-5p, miR-152 and miR-653 were involved. Analysis of the pathways indicated that only miR-146b-5p and miR-152 were directly related to the pancreatic cancer pathway (hsa05212: pancreatic cancer). Finally, 5 prognosis-related KEGG pathways were obtained after enrichment analysis ( Table 2), 3 of which overlapped with miR-146b-5p-and miR-152-related pathways: hsa04350: TGF-beta signaling pathway, hsa04360: axon guidance and hsa04530: tight junction. A ceRNA network of the above 3 pathways was constructed, consisting of 4 pathway nodes, 2 upregulated miRNAs, 5 downregulated circRNAs, 37 downregulated mRNAs and 3 upregulated mRNAs (Fig. 4. C). The downregulation of hsa_circ_0006502 was veri ed in pancreatic cancer cell lines The expression levels of the 8 most differentially expressed circRNAs, as mentioned above, were tested in 5 pancreatic cancer cell lines. Hsa_circ_0006502, hsa_circ_0000919, hsa_circ_0079385, hsa_circ_0092367, hsa_circ_0002191 and hsa_circ_0013587 were expected to be downregulated, while hsa_circ_0000376 and hsa_circ_0008240 were expected to be overexpressed in cancer cells compared to normal pancreas cells. The 8 circRNAs did not behave as expected in all ve cancer cell lines, except hsa_circ_0006502, which was con rmed to be decreased in all ve pancreatic cancer cell lines, as predicted. The qRT-PCR results for circRNA veri cation were acceptable, considering the heterogeneity of cancer cells (Fig. 5).

Hsa_circ_0006502 Functions As A Tumor Suppressor In Pancreatic Cancers
Stable overexpression of hsa_circ_0006502 was established in the CFPAC-1 and AsPC-1 cell lines. The level of hsa_circ_0006502 was markedly increased in CFPAC-1 and AsPC-1 cells transfected with plasmid pLCDH_ciR_hsa_circ_0006502 (supplementary Fig. 3). Compared with that in the control group, CFPAC-1 and AsPC-1 cell proliferation was signi cantly suppressed in the hsa_circ_0006502 overexpression group (Fig. 6. A). The apoptosis results also revealed that CFPAC-1 and AsPC-1 apoptotic cells were markedly increased by hsa_circ_0006502 overexpression (Fig. 6. B). Transwell assays demonstrated that hsa_circ_0006502 overexpression markedly inhibited the migratory and invasive capabilities of CFPAC-1 and AsPC-1 cells (Fig. 6. C). Overall, hsa_circ_0006502 acted as a tumor suppressor in CFPAC-1 and AsPC-1 cells.

Discussion
By integrating different types of RNAs (circRNA, miRNA, and mRNA) from multiple public databases, RNAs with signi cantly altered expression levels in tumor tissues compared to normal tissues were identi ed, based on which a ceRNA regulatory network was constructed. In total, 798 differentially expressed genes between pancreatic cancer tissues and adjacent tissues were obtained, consisting of 85 circRNAs, 19 miRNAs and 694 mRNAs, from all 14 databases. The results from the GO analysis indicated that the differentially expressed genes were mostly involved in biological processes such as apoptosis, programmed cell death, tube development and growth regulation. KEGG pathway enrichment analysis suggested that these genes were enriched in pathways such as tight junction, starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism, MAPK signaling and Notch signaling. These biological processes and enriched pathways have crucial functions in many different solid tumors and provide novel and detailed insights into the molecular mechanism of pancreatic cancer that could help develop therapeutic strategies in pancreatic cancers.
Three miRNAs and 75 mRNAs signi cantly related to pancreatic cancer prognosis were identi ed, and a prognosis-related ceRNA network was constructed. MiR-146b-5p has important and complex functions and its dysregulation been discovered in a variety of human malignancies. It is overexpressed in some malignant tumors, such as colorectal cancer [23], thyroid cancer [24], renal cell carcinoma [25], and osteosarcoma [26], and acts as an onco-miRNA, promoting cancer growth and metastasis. In other cancers, such as non-small cell lung cancer [27], gallbladder carcinoma [28], and glioma [29], miR-146b-5p functions as a tumor suppressor, inhibiting cancer development. Although miR-146b-5p was reported to help suppress the proliferation of pancreatic cancer cells [30], miR-146b-5p expression was inconsistent in pancreatic cancer development, as it was downregulated in GSE25820 but upregulated in GSE60978, and in TCGA_PAAD, it was overexpressed in high-grade cancers compared with low-grade cancers [31]. In our present study, miR-146b-5p was revealed to be overexpressed in pancreatic cancers by meta-analysis, and its lower expression usually indicated a longer survival time. Similar to miR-146b-5p, miR-653 plays different roles in different cancers. It appears to function as a tumor suppressor in most malignant tumors, such as renal cell carcinoma, melanoma, gastric cancer, cervical cancer, breast cancer, Wilms' tumor and non-small cell lung cancer [32][33][34][35][36][37][38], but it was also reported to act as an onco-miRNA in prostate cancer and bladder cancer [39,40]. Regardless, our meta-analysis indicated that miR-653 was downregulated in pancreatic cancer, suggesting that it might suppress pancreatic cancer. MiR-152 is frequently downregulated and has a tumor-suppressive function in different solid tumors [41], but it was predicted to be upregulated in our present study, similar to other datasets [31]. Such inconsistency re ects the limitations of bioinformatics analyses and illustrates the importance of experimental veri cation.
Bioinformatics analyses can help to identify potentially valuable questions, but more in-depth investigation is needed to answer such questions.
CircRNAs have a stable circular structure and long half-lives, which makes them suitable for noninvasive detection and useful as biomarkers and therapeutic targets. Emerging studies have demonstrated that circRNAs can be signi cantly differentially expressed between pancreatic cancer tissues and adjacent tissues and are involved in cancer pathogenesis with promising clinical value in diagnosis and treatment [42][43][44][45]. Here, the 8 circRNAs included in the prognosis-related regulatory network were veri ed in ve pancreatic cancer cell lines. The expression levels of some of the circRNAs did not behave as predicted by the bioinformatics analysis across all cancer cell lines. Interacting closely with miR-146b-5p, hsa_circ_0006502 was con rmed to be downregulated in all pancreatic cancer cell lines, consistent with the prediction by the bioinformatics analysis. Functional tests demonstrated that hsa_circ_0006502 inhibited cancer cell growth, migration and invasion and acted as a tumor suppressor. However, further studies are needed to reveal the underlying mechanism of tumorigenesis and the clinical value of hsa_circ_0006502 and its associated regulatory network.
In conclusion, a key prognosis-related regulatory network in pancreatic cancers was constructed through bioinformatics analysis of public RNA databases such as NCBI GEO and EBI ArrayExpress. CircRNAs with signi cantly aberrant expression levels included in the network were veri ed by qRT-PCR in ve pancreatic cancer cell lines. Hsa_circ_0006502 was identi ed as a new tumor suppressor that inhibits cancer cell growth, migration and invasion and might serve as a potential therapeutic agent.

Declarations
Ethics approval and consent to participate Not applicable