Network analyses of circular RNAs expression profiles and their hub genes in pancreatic ductal adenocarcinoma

Pancreatic ductal adenocarcinoma (PDAC) is one of the most common malignant tumor in digestive system. CircRNAs involve in lots of biological processes through interacting with miRNAs and their targeted mRNA. We obtained the circRNA gene expression profiles from Gene Expression Omnibus (GEO) and identified differentially expressed genes (DEGs) between PDAC samples and paracancerous tissues. Bioinformatics analyses, including GO analysis, KEGG pathway analysis and PPI network analysis, were conducted for further investigation. We also constructed circRNA‑microRNA-mRNA co-expression network. A total 291 differentially expressed circRNAs were screened out. The GO enrichment analysis revealed that up-regulated DEGs were mainly involved metabolic process, biological regulation, and gene expression, and down-regulated DEGs were involved in cell communication, single-organism process, and signal transduction. The KEGG pathway analysis, the upregulated circRNAs were enriched cGMP-PKG signaling pathway, and HTLV-I infection, while the downregulated circRNAs were enriched in protein processing in endoplasmic reticulum, insulin signaling pathway, regulation of actin cytoskeleton, etc. Four genes were identified from PPI network as both hub genes and module genes, and their circRNA‑miRNA-mRNA regulatory network also be constructed. Our study indicated possible involvement of dysregulated circRNAs in the development of PDAC and promoted our understanding of the underlying molecular mechanisms.


Introduction
Pancreatic ductal adenocarcinoma (PDAC), which is the most common form of pancreatic carcinoma, remains the fourth highest mortality rate malignant neoplasms in the United States. More specifically, its estimated deaths will be 43,090 in 2017 [1].
The prognosis of PDAC is poor. It is a highly metastatic tumor and therapy resistant. Its survival rate is approximately 25% of 1-year and 5% of 5-year [2]. Beside the dietary habit and environmental factors, PDAC is also related with genetic basis. The family history of pancreatic carcinoma is predisposing to increased risk in direct relatives [3].
The early diagnostic is very hard due to lack of characteristic symptoms and cancer late stage during first oncological visit. It is usually not adequate to treat this disease by surgical resection [4]. There are still need potent specific biomarkers for early diagnosis. System analyses of PDAC may provide better understanding to this disease and contribute to developing more effective approaches. Circular RNAs (circRNAs) was first declared in 1976 [5]. However, circRNAs were almost completely neglected since then, because it is difficult to be detected using traditional molecular techniques, and hard to be mapped into genomes [6]. With the development of sequencing technologies and bioinformatics methodologies, circRNAs have become the focal research in the biological field [7].
CircRNAs are unusually regarded as stable molecules. Their 3' and 5' ends are linked together covalently [8], and thus lack of a free end, which prevents them depredated as the conventional RNA.
The special features of circRNA, such as stability, abundance, and tissue-specific expression pattern, make them very attractive for clinical research, as they can be employ as new ways in diagnosis [4].
CircRNAs can repress the function of microRNA (miRNA) by binding to miRNA as sponges [9]. They involve in lots of biological processes by interacting with miRNAs and their target mRNA. CircRNAs have been reported in several types of diseases, including Alzheimer's disease [10], colorectal cancer [11], atherosclerotic vascular disease [12], and osteoarthritis [13].
To explore the underlying molecular regulation mechanism of circRNAs in PDAC, circRNAs microarray was used to detect the differential expressed circRNAs between PDAC tissues and paracancerous tissues. Our results suggested that the dysregulated expression of circRNAs may play a role in transformation of PDAC.

Microarray data
The circRNA expression profiles were obtained from Gene Expression Omnibus (GEO), where the accession number is GSE79634. It contains 40 samples, including twenty fresh-frozen PDAC samples and twenty paracancerous tissues. The circRNA expression levels of these samples were measured using Agilent-069978 Arraystar Human circRNA microarray V1. Total RNA was extracted from collected samples. After labeling and hybridizing, the acquired array images were analyzed using Agilent Feature Extraction software (version 11.0.1.1). Quantile normalization and subsequent data processing were performed using the R software package.

Identification of differentially expressed genes
Differentially expressed genes (DEGs), which were significantly differentially expressed among the PDAC samples and the paracancerous tissues, were identified by using the LIMMA package in R language. LIMMA provides several p-value adjustment options to correct for the occurrence of false positive results. It was also known as multiple-testing corrections, where the Benjamini & Hochberg false discovery rate method was selected by default in the package. We defined adjusted p-values ≤0.05 and absolute value of fold change ≥ 2 as the cutoff to determine the statistically significance of gene expression difference.

GO enrichment analysis
Gene Ontology (GO) is a common useful framework for the model of biology. The GO project provides a computational representation of our evolving knowledge on how genes encode biological functions.
It has defined formal ontologies that represent over 40,000 biological concepts at the molecular, cellular and system levels, to describe gene function and relationships between them.

KEGG pathway analysis
Kyoto encyclopedia of genes and genomes (KEGG) database resource is an encyclopedia of genes and genomes for understanding high-level functions and utilities of the biological system from molecularlevel information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies. It assigns functional meanings to genes and genomes both at the molecular and higher levels. Molecular level functions are maintained in the KO (KEGG Orthology) database, while higher-level functions are defined in the forms of KEGG modules, KEGG pathway maps and BRITE hierarchies [14].

PPI network analysis
It is widely accepted that the interactions and partnerships between proteins are play important roles in biological regulation. The protein-protein interaction networks, which are constructed from either known or predicted interactions, provide us with topological summary of our research subjects.
A larger number of online resources contain this protein-protein interaction information. These interactions were generated either from experimental techniques [15,16] or from computational prediction [17][18][19]. A group of databases integrate both known and predicted interactions, to provide much high comprehensiveness and coverage. Typical databases include GeneMANIA [20], FunCoup [17], I2D [21], and STRING [22].
Highly complex biological functions are often undertaken by tightly interacted groups of proteins. PPI network was constructed in Cytoscape, and the significant modules were identified using the plug-in Molecular Complex Detection (MCODE).

Target prediction and co-expression network
To further investigate the functional roles of circRNA, the miRNAs regulated by these circRNA and corresponding target genes were extracted from CIRCNET database. The circRNA-microRNA-mRNA coexpression networks were constructed based on this information by using Cytoscape.

Differentially expressed genes
Differentially expressed genes (DEG) were identified from microarray data, using the criteria of adj. P.

Enrichment analysis of DEGs
The Gene Ontology (GO) enrichment analyses were performed on these 128 upregulated and 163 downregulated circRNAs respectively. Table 2 was the summary of GO enrichment analysis results, which listed the number of enriched GO terms at different significant level (P<0.10 or P<0.05). Table   3 and Table 4 were top 10 GO terms (P<0.05), respectively. GO enrichment analysis revealed that numerous genes were involved in the biological processes. For example, according to the gene ontology enrichment analysis results, the up-regulated circRNAs were enriched (P <0.05) in 80 biological processes, including biological regulation, metabolic process, and gene expression, etc.
Meanwhile, the downregulated circRNAs were enriched (P <0.05) in 151 biological processes, including single-organism process, cell communication, and signal transduction, etc. These processes were associated with PDAC.

KEGG pathway analysis
KEGG pathway analysis results were listed in Table 5. The cutoff used for this analysis was P<0.05.
According to these results, the upregulated circRNAs were enriched in 2 pathways including cGMP-PKG signaling pathway, and HTLV-I infection, while the downregulated circRNAs were enriched in 10 pathways, including Protein processing in endoplasmic reticulum, insulin signaling pathway, regulation of actin cytoskeleton, etc.

PPI network analysis and hub genes screening
The PPI (protein-protein interaction) network of the dysregulation genes were constructed by Cytoscape software, where the interaction information was extracted from the STRING database. As shown in Figure 2A, the red triangle nodes represent upregulated genes and green circular nodes represent downregulated genes, respectively.
The PPI network was further analyzed by screening its hub genes and identifying significant modules.
In the PPI networks, the nodes with high degree are defined as hub genes. The top 10 hub genes were ACACA, ACACB, GSK3B, LRRK2, PRKG1, UBQLN1, FOXO1, AXIN1, MTHFD1L and PTK2. The most significant hub proteins in the PPI network were ACACA (degree = 19) and ACACB (degree = 18). We also analyzed the PPI network using MCODE, which is a Cytoscape plug-in that finds clusters (or modules, highly interconnected regions) in a network. In PPI network, clusters are often protein complexes and parts of pathways. There were 3 modules identified in this PPI network ( Figure 2A).
3.5. CircRNA-microRNA-mRNA co-expression network As shown in Table 6, there were four genes identified as both hub genes and module genes in the PPI network analysis, including ACACA, MTHFD1L and LRRK2. The circRNA-microRNA-mRNA co-expression networks were constructed based on these genes.
The miRNAs regulated by these circRNA and corresponding target genes were extracted from CIRCNET database. Due to limited space, no more than five top significant regulated miRNAs of each circRNA were reserved for the construction of co-expression network, and so do their target genes. In the circRNA-miRNA-mRNA regulatory network, there also existed closed loops. That is, circRNA regulated miRNA through its specific transcript, and miRNA may also regulate circRNA in turn, thus forming a closed loop structure. For example, in Figure 4, ACACB regulated miR-6849-3p through circ-ACACA.15, and miR-6849-3p in turn inhibits the expression of ACACB, thereby forming a closed loop.

Discussion
Pancreatic ductal adenocarcinoma is an aggressive malignant neoplasm with extremely high mortality rate, mainly because of difficult early diagnostic and poor prognosis. Current available therapeutic methods are almost ineffective [23]. The potential improvement of the diagnosis may exist in the strictly molecular biomarkers [24]. CircRNA is a new special kind of RNAs. Recent evidences revealed that circRNAs can function as miRNA sponges and affect disease by regulating gene expression.
In the present study, we employed bioinformatics technology to investigate the underlying molecular mechanism of PDAC. Differentially expressed profiles of circRNAs in PDAC tissues were observed and validated compared with the paracancerous tissues, indicating possible involvement of these dysregulated circRNAs in the development of PDAC.
By using microarray data, we observed a total of 291 circRNAs differentially expressed between PDAC tissues and paracancerous tissues, including 128 up-regulated circRNAs and 163 down-regulated circRNAs. Among these differentially expressed genes, hsa_circ_0006220, hsa_circ_0000977 and hsa_circ_0001666 were upregulated with top magnitudes, while hsa_circ_0000518, hsa_circ_0072088 and hsa_circ_0005273 were downregulated with top magnitudes. The differential expression levels of these genes may be related to the involvement in the transcription level regulation on PDAC.
Cumulative evidence has demonstrated that co-expression genes usually participate in similar biological process. In order to better understand the interactions of these differentially expressed genes, we further applied GO and KEGG pathway analysis on them.
The GO enrichment analysis results revealed that up-regulated DEGs were mainly involved metabolic process, biological regulation, and gene expression, and down-regulated DEGs were involved in cell communication, single-organism process, and signal transduction.
As is well-known, the symptoms of the patients with cancer, such as cachexia, wasting syndrome, is related to metabolic disorders. Recent studies showed that, in the case of caloric deficiency, tumor induced changes in host metabolism, leading to tumor immunosuppression. For example, in precachexia mice with PDA, IL6 reduced liver ketogenic ability by inhibiting PPAR alpha. When these mice challenged with caloric deficiency, it would lead to hypoketonemia, and induce the elevation of glucocorticoid levels. A variety of tumor immune pathways was inhibited by above stress response, and thus promoted the progression of pancreatic cancer [25]. Furthermore, poor survival rates in patients with PDAC were attributed to abnormally high levels of metastasis. Study showed that the metabolic reprogramming of tumor infiltrating macrophages played a key role in PDAC metastasis [26]. In addition, autophagy is required to maintain energy homeostasis by degrading unnecessary cellular components and molecules. In recent years, it has been found that autophagy can regulate the metabolism of tumor. Autophagy plays a crucial role in the maintenance of tumor survival. PDAC cells require autophagy and glutamine transporters to maintain intracellular glutamine levels [27].
Therefore, we hypothesized that the high expressed circRNAs may promote the metastasis of pancreatic cancer through various metabolic process, such as, stress response, metabolic reprogramming and glutamine metabolism.
PDAC is a highly metastatic tumor with poor prognosis, and intercellular communication is the key to metastasis [28]. Research showed that PDAC-derived exosomes induced liver pre-metastasis niche formation in nude mice, and increased the load of liver metastasis. Kupffer cells uptaked PDACderived exosomes, resulting in transforming growth factor β (TGF-β) and upregulation of fibronectin production by hepatic stellate cells. This fibrotic microenvironment increases the number of bone marrow-derived macrophages. The authors found that macrophage migration inhibitory factor (MIF) was highly expressed in PDAC derived exosomes, and MIF promote metastasis, and liver metastasis before niche formation. It can be seen that the MIF of the exocrine body is initially transferred to the liver, which may be a prognostic indicator of PDAC liver metastasis [28]. used the KEGG analysis and found that many endoplasmic reticulum (ER)-associated proteins were altered in solid pseudopapillary tumor of the pancreas (SPTP), suggesting that endoplasmic reticulum stress may play an important role in SPTP tumorigenesis. They also verified 7 proteins in the pathway by immunohistochemistry, including ERO1LB, TRIM1, BIP, SEC61B, P4HB, GRP94 and, PDIA4. Six of these proteins (except SEC61B) were consistent with the LC-MS/MS results [29].
Our KEGG analysis results also showed that the down regulated circRNAs of pancreatic cancer were enriched in insulin signaling pathway. It has been reported that hyperinsulinemia was a clear risk factor for pancreatic cancer, and an increase in hyperinsulinemia associated with obesity and type 2 diabetes heralded an increased incidence of cancer. The role of insulin in pancreatic cancer progression was worth to be further studied [30]. They suggested that the overexpression of insulin signaling contributes to the survival and proliferation of human immortalized pancreatic ductal cells and metastatic pancreatic cancer cells, but does not affect normal pancreatic ductal cells. They also suggested that cell survival associated insulin pathway may be involved in the progression of pancreatic cancer [30]. There were also literatures about insulin signaling pathway and chemotherapeutic resistance of pancreatic cancer. Ireland et al. found that tumor associated macrophages (TAM) and myofibroblasts activated insulin /IGF receptor by secreting IGF1-1 and 2, and directly promote chemoresistance of pancreatic cancer cells resistance. Analysis of biopsies from patients with pancreatic cancer showed that 72% of the patients' tumor cells had expressed activated insulin / IGF receptor. Therefore, inhibition of IGF might contribute to the treatment of PDAC patients with chemotherapy resistance, and activation of insulin /IGF1R might be used as a biomarker for the identification of chemotherapy resistance [31]. In addition, adaptor protein CrkII played an important role in the proliferation and invasion of some malignant tumors. However, the specific mechanism of insulin like growth factor 1 (IGF-1) -CrkII signal induced proliferation of PDAC was not clear. Studies have further demonstrated that the CrkII in insulin signaling pathway mediated IGF-1 signaling and affected PDAC through Erk1/2 and Akt pathway, suggesting that CrkII gene and protein could be used as an effective therapeutic targets for PDAC [32].
We further applied PPI network analysis to these aberrant expression genes. Ten hub genes with the top degree were obtained from the PPI network. We also identified three top significant modules by clustering. A total of eleven module genes were contain in these three module. Four genes were existed both as hub genes and module genes, including ACACA, MTHFD1L and LRRK2.
We predicted the target genes of these four circRNA genes, and constructed circRNA-miRNA-Gene regulatory networks for them. On one hand, These genes regulated the expression of miRNAs by transcribing a variety of circRNAs, and further regulated the target genes through miRNAs. According to these regulatory networks, the more circRNAs, the more complex the regulation of miRNA and target genes. There existed multiple miRNAs regulated by the same circRNA. There also existed the same miRNA regulated by several circRNAs of a specific gene. In addition, circRNAs of different genes might regulate the same miRNA. On the other hand, On the other hand, miRNA also could bind to the target gene, and then regulated circRNAs via negative feedback, thus forming a closed loop. For example, circ-ACACB.74, circ-ACACB.25, and circ-MTHFD1L.101 targeted miR-508-5p simultaneously, and miR-508-5p regulated the expression of ACACB and MTHFD1L through the negative feedback.
Therefore, bioinformatic analyses showed that circRNA targeting miR-508-5p might involve in the extremely complex regulation mechanism of PDAC. There exist literatures that miR-508-5p involves in the regulation mechanisms of liver gastric cancer and hepatocellular carcinoma. Shang Y. et al.
reported that miR-508-5p play an important role in multidrug resistance (MDR) of gastric cancer by targeting ABCB1 and ZNRD1 [33]. Further research suggested that the underlying mechanism was related to the miR-27b/CCNG1/P53/miR-508-5p axis, and it was possible to reverse MDR of gastric cancer by restoring the expression levels miR-27b and miR-508-5p [34]. Yan H. et al. found out that miR-508-5p was aberrant expressed in hepatocellular carcinoma (HCC) [35]. Wu et al. suggested that miR-508-5p might act as a tumor suppressor of HCC progression by targeting MESDC1 [36]. Similar to gastric cancer and hepatocellular carcinoma, pancreatic adenocarcinoma is also a disease belong to the digestive system. We have reason to believe that miR-508-5 p may also be in pancreatic adenocarcinoma plays an important role. It is worth to study the problem of whether miR-508-5p participates in PDAC.

Data Availability
The microarray data used to support the finding of this study have been deposited in the GEO expression associated with paclitaxel-induced apoptosis in hepatocellular carcinoma cells.