Potential ceRNA networks involved in autophagy suppression of pancreatic cancer caused by chloroquine diphosphate: A study based on differentially‑expressed circRNAs, lncRNAs, miRNAs and mRNAs.

Autophagy has been reported to be involved in the occurrence and development of pancreatic cancer. However, the mechanism of autophagy‑associated non‑coding RNAs (ncRNAs) in pancreatic cancer remains largely unknown. In the present study, microarrays were used to detect differential expression of mRNAs, microRNAs (miRNAs), long ncRNAs (lncRNAs) and circular RNAs (circRNAs) post autophagy suppression by chloroquine diphosphate in PANC‑1 cells. Collectively, 3,966 mRNAs, 3,184 lncRNAs and 9,420 circRNAs were differentially expressed. Additionally, only two miRNAs (hsa‑miR‑663a‑5p and hsa‑miR‑154‑3p) were underexpressed in the PANC‑1 cells in the autophagy‑suppression group. Furthermore, miR‑663a‑5p with 9 circRNAs, 8 lncRNAs and 46 genes could form a prospective ceRNA network associated with autophagy in pancreatic cancer cells. In addition, another ceRNA network containing miR‑154‑3p, 5 circRNAs, 2 lncRNAs and 11 genes was also constructed. The potential multiple ceRNA, miRNA and mRNA associations may serve pivotal roles in the autophagy of pancreatic cancer cells, which lays the theoretical foundation for subsequent investigations on pancreatic cancer.

Introduction

Pancreatic cancer is a highly aggressive and fatal malignancy, with stably high incidence and mortality for the past 30 years globally (1-4). Patients with pancreatic adenocarcinoma account for ~85% of all patients with pancreatic cancer and have a poor prognosis, compared with patients with other common solid tumors (2,5-9). Pancreatic cancer is primarily treated through surgical resection, whereas postoperative patients and patients with advanced disease are frequently subjected to adjuvant therapies with radiotherapy, chemotherapy, and a combination of radiotherapy and chemotherapy that are associated with high recurrence rates and poor treatment effectiveness (10). The recent development of and research into molecular targeted drugs have great significance for the early diagnosis and prognosis of malignant tumors, and have ushered a new era of tumor therapy (11-15). However, the overall treatment effect of these drugs remains poor, and further investigation is required (16).

Autophagy is a process by which eukaryotic cells undergo self-degradation of intracellular damaged macromolecular proteins and organelles (17,18). Autophagy at the physiological level is necessary to maintain the stability of the internal environment, and the occurrence of a variety of diseases, including inflammation, tumors and degenerative diseases, is frequently accompanied by an abnormal autophagic level (19-25). Currently, well-characterized, autophagy-associated proteins include Beclin 1, microtubule-associated-protein-1-light-chain-3 (LC3) and P62, and well-studied pathways include mechanistic target of rapamycin (mTOR) and phosphoinositide 3-kinase (PI3K) pathways (26,27).

Recently, the important role of autophagy in tumors has been recognized (19,28). Studies demonstrated that pancreatic cancer cells have an increased autophagy level, compared with other tumor cells (29-34). The growth of PDAC had a distinct dependence on autophagy in vivo and in vitro (35). Inhibition of autophagy is effective for the treatment of pancreatic cancer through animal experiments (36). Another study revealed that autophagy can be activated by gemcitabine and ionizing radiation in the treatment of pancreatic cancer cells, and activated autophagy serves a role in cancer suppression (37). However, whether autophagy promotes the occurrence, development and prognosis of pancreatic cancer is controversial (38). Therefore, the role of autophagy in the occurrence, development, prognosis and treatment of pancreatic cancer and its mechanism requires further examination.

mRNA is a single-stranded RNA that is transcribed using one of the DNA chains as the template and carries the genetic information that guides protein synthesis. Thus, mRNA occupies an important position in the heredity dogma (39). mRNA serves a key role in various human diseases and is involved in autophagy pathways (40,41); for example, signal transducer and activator of transcription 3 could inhibit autophagy in pancreatic cancer cells (42) and ubiquitin specific peptidase 1 regulates autophagy by targeting Unc-51 like autophagy activating kinase 1 (ULK1) (43). Non-coding RNA (ncRNA) refers to RNA that does not encode a protein and includes ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), long ncRNA (lncRNA), microRNA (miRNA) and circular RNA (circRNA) (44-50). The involvement of lncRNAs and miRNAs in cancer has been studied intensively, whereas the potential significance of circRNAs in tumors has attracted attention over the last two years (51). The post-transcriptional levels of circRNAs serve a notable regulatory role in gene expression (52). Furthermore, the sponge effect of circRNAs on miRNAs is an important mechanism (53). circRNAs can contain multiple miRNA binding sites, and attract miRNAs effectively to reduce miRNA-mediated mRNA inhibition (54-56). Based on the interaction between these three types of ncRNAs and the current research status in pancreatic cancer, it was considered that an investigation of the effect of lncRNAs and circRNAs as competitive endogenous RNAs (ceRNAs) on the targeting of miRNAs in the development of pancreatic cancer is beneficial.

Recently, the number of autophagy-associated lncRNAs has been determined to be limited (57-60), and investigations on circRNAs and their associations with autophagy have not been reported to date. However, due to autophagy being involved in the occurrence and development of a variety of diseases, including breast (61,62), gastric (63) and lung cancer (64), in the present study, the human pancreatic cancer cell line PANC-1 and gene chip technology was used to detect differential expression of mRNAs, miRNAs, lncRNAs and circRNAs under different autophagy levels, and to investigate the genes associated with autophagy in pancreatic cancer and their underlying molecular mechanisms.

Materials and methods

Cell treatment

The human pancreatic cancer cell line PANC-1 was obtained from the Cell Laboratory of Chinese Academy of Sciences (Shanghai, China) and cultured in complete Dulbecco’s modified Eagle’s medium (DMEM; Gibco; Thermo Fisher Scientific, Inc., Waltham, MA, USA) at 37°C. Due to chloroquine diphosphate being frequently used as a classic autophagic inhibitor (65,66), the following concentration gradient of chloroquine diphosphate (Sigma-Aldrich; Merck KGaA, Darmstadt, Germany) was designed as recommended by the manufacturer’s protocols: 0, 25, 50, 75 and 100 µM. When the number of cells in the culture bottle grew to ~5×106, the original medium was replaced by a complete DMEM containing chloroquine diphosphate. Following treating the cells for 12 h at 37°C with the drug, the cellular proteins were immediately extracted using Radioimmunoprecipitation Assay Lysis Buffer (Beyotime Institute of Biotechnology, Haimen, China).

Western blotting

The lysis solution (radioimmunoprecipitation assay buffer:phenylmethylsulfonyl fluoride; 100:1, Beyotime Institute of Biotechnology) was prepared. The protein concentration in the lysate was determined with a bicinchoninic acid kit (Beyotime Institute of Biotechnology), according to the manufacturer’s instructions. Western blotting was performed according to routine procedures, and conducted with a 12% separating gel and a 5% stacking gel. A total of 20 µg protein was added to the glue holes. The electrophoresis condition was set to 120 Ma for 1 h. Following the electrophoresis, the protein was transferred to the polyvinylidene (PVDF) membrane (Beyotime Institute of Biotechnology). Subsequently, the PVDF membrane was placed into the Bull Serum Albumin blocking buffer (Beyotime Institute of Biotechnology) for shaking at 37°C for 2 h. Primary antibodies were incubated at 4°C overnight and secondary antibodies were incubated at room temperature for 2 h. The primary antibodies were LC3 β-specific rabbit polyclonal (1:1,000; cat. no. 18725-1-AP) and P62 mouse monoclonal antibodies (1:2,000; cat. no. 66184-1-Ig), which represented the state of autophagy, and the GAPDH horseradish peroxidase-labeled internal reference antibody (1:10,000; cat. no. HRP-60004), andhe secondary antibodies used were goat anti-mouse IgG (H+L), allopregnanolone conjugate (1:1,000; cat. no. SA00002-1) (all from ProteinTech Group, Inc., Chicago, IL, USA). The gray scale values of the protein bands from the raw image were determined using the ImageJ 1.48 (National Institutes of Health, Bethesda, MD, USA) and the formula Final gray scale value (G) = [G × (target band) − G × (background)] / [G × (internal reference band of the group) − G × (background)]. The data were analyzed using SPSS software 22 (SPSS, Inc., Chicago, IL, USA), and the results were plotted using the GraphPad Prism 5 software (GraphPad Software, Inc., La Jolla, CA, USA).

Sample grouping

Based on the results of the western blotting experiment, the human pancreatic cancer PANC-1 cells treated with 100 µM chloroquine diphosphate were selected as the autophagic inhibition group, and the cells cultured in normal DMEM medium were selected as the control group. Subsequent experiments were compared between the autophagic inhibition group and the control group. Following chloroquine diphosphate treatment, total RNA was prepared using the TRIzol® reagent (Beyotime Institute of Biotechnology) with three replicates per group.

ceRNA microarray and detection, and statistical analysis of the miRNA chip results

The quality control for the gene analyses in the present study was performed by Shanghai Biotechnology Corporation (Shanghai, China). The starting sample of the microarray test was total RNA, which was analyzed with the NanoDrop ND-2000 spectrophotometer and the Agilent Bioanalyzer 2100 (Agilent Technologies, Inc., Santa Clara, CA, USA) to conduct the quality control assay. Only RNA samples that passed the quality control proceeded to the subsequent microarray experiments. The qualifying sample standards were that the RNA integrity number of each sample was ≥7.0 and that the 28S/18S was ≥0.7.

The ceRNA microarray detection assay included the detection of three types of RNA (lncRNA, circRNA and mRNA). The cutoff values of the differentials were all set to a fold change (FC)>2 or FC<0.5 and P<0.05. The assays and data analyses are described below. The analyses primarily included the normalization of raw data, sample association analysis, screening of genes with differential expression, Gene Ontology (GO) enrichment analysis of the differentially-expressed mRNAs, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, PANTHER pathway analysis of the differentially-expressed mRNAs, prediction of lncRNA and circRNA target genes, and prediction of lncRNA and circRNA adsorption of miRNAs. GO and KEGG pathway were generated by DAVID (https://david.ncifcrf.gov/) while PANTHER was analyzed using the WebGestalt 2017 tool (http://www.webgestalt.org/option.php). The clinical roles of genes were analyzed through the online database of Gene Expression Profiling Interactive Analysis (GEPIA; http://gepia.cancer-pku.cn/) Which is based on the data obtained from The Cancer Genome Atlas and GTEx. The miRNA microarray analyses primarily included the normalization of raw data, sample association analysis, screening of differentially-expressed miRNAs and prediction of miRNA target genes, which was predicted by 12 online miRNA target prediction databases [DIANA-microTv4.0 (67), DIANA-microT-CDS (68), miRanda-rel2010 (69), mirbridge (70), miRDB4.0 (71), miRMap (72), miRNAMap (73), PicTar2 (74), PITA (75), RNA22v2 (76), RNAhybrid2.1 (77) and Targetscan6.2 (78)] (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/). The crossed genes from target prediction and differentially-expressed mRNAs were depicted in Compendia expression profiles (Novartis; http://software.broadinstitute.org/gsea/msigdb/annotate.jsp). Network visualization was performed in Cytoscape 3.5.1 (The Cytoscape Consoritum, New York, NY, USA).

The raw data of ceRNA microarray has been deposited in the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/; accession number GSE115517).

Tissue samples

A total of 31 formalin-fixed paraffin-embedded (FFPE) tissues, including 18 cases of Pancreatic ductal adenocarcinoma and 13 paracancerous tissues, were collected from the Department of Pathology of the First Affiliated Hospital of Guangxi Medical University (Nanning, China) between June 2015 and June 2018. The present study was approved by the Ethics Committee of the First Affiliated Hospital of Guangxi Medical University. All cases were from patients who have not been treated with chemotherapy or radiation prior to resection and have signed informed consent. There were 11 male and 7 female aged 29-77 years (mean age, 56 years) in 18 patients with adenocarcinoma.

RNA extraction and Real-time fluorescent quantitative PCR (RT-qPCR)

According to the user guides provided by the manufacturer, the total RNA of tissues was extracted using an E.Z.N.A.® FFPE RNA kit (Omega Bio-Tek, Inc., Norcross, GA, USA). The complementary DNA (cDNA) was reverse transcribed through the kit of miRNA First Strand cDNA Synthesis [Tailing Reaction (79,80)] (Sangon Biotech Co., Ltd., Shanghai, China), and RT-qPCR was conducted by applying a MicroRNAs qPCR kit [SYBR® Green method (81,82)] (Sangon Biotech Co., Ltd.), according to the manufacturer’s protocols, on an ABI Prism 7500 (Applied Biosystems; Thermo Fisher Scientific, Inc.). The temperature of pre-denaturation and denaturation were both set at 95°C, and the temperature of annealing/extension was set at 60°C. Pre-denaturation was conducted for 10min, denaturation for 15 sec and annealing/extension for 60 sec. The number of cycles was set to 40. Subsequently, the expression of miRNA relative to U6 was calculated via the 2−ΔΔCq method on identical samples (83). The PCR primers included: i) miR-663a-5p for ward, 5′-ATAGGCGGGGCGCCGCGGGAC-3′; ii) miR-154-3p forward, 5′-CCGGGAATCATACACGGTTGA CCTATT-3′; and iii) U6 forward and universal PCR reverse primer were attached to the kits. Their primer sequences are 5′-CTCGCTTCGGCAGCACA-3′ and 5′-AACGCTTCACGA ATTTGCGT-3′, respectively.

Retrieval of data from GEO

Within GEO, ‘pancreas’ or ‘pancreatic’, and ‘adenocarcinoma’, ‘carcinoma’, ‘cancer’, ‘neoplasm’, ‘tumor’, ‘tumour’, ‘neoplas*’, ‘malignan*’, ‘PDAC’, ‘OR’, ‘PAAD’ or ‘PC’ were employed as a search strategy in GEO to determine the expression of miR-663a-5p and miR-154-3p. STATA 12.0 (StataCorp LLC, College Station, TX, USA) was applied to estimate the pooled effects, heterogeneity and publication bias.

Statistical analysis

The data of GEO and RT-qPCR were calculated using SPSS 22 software for mean and standard deviation (SD). Continuous data are presented as mean ± SD. The forest plots were produced by Stata 12.0. The heterogeneity test was used to analyze the existence of heterogeneity and the source of the occurrence of heterogeneity (I2>50% or P<0.05 is the existence of heterogeneity). Funnel plots were used for the analysis of publication bias. The comparison between the two groups was performed using two-sample Student’s t-test. The method for differential gene expression analysis in GEPIA was one-way analysis of variance, using pathological stage as variable for calculating differential expression. GEPIA performs overall survival analysis based on gene expression. GO analysis was performed in DAVID tool which used Fisher’s exact test and the techniques of Kappa statistics (84). PANTHER analysis was performed in WebGestalt, which used Hypergeometric and Fisher’s exact tests. P<0.05 was considered to indicate a statistically significant difference.

Results

Autophagy following chloroquine diphosphate treatment

The present study design of the whole investigation is depicted in Fig. 1. The LC3-II expression level was significantly reduced in the control group, compared with the group treated with 25 or 100 µM chloroquine diphosphate (P<0.05). However, the P62 expression level was significantly increased in the group treated with 100 µM chloroquine diphosphate, compared with the control group or the groups treated with 25, 50, 75 or 100 µM chloroquine diphosphate (P<0.05). Thus, the 100 µM concentration of chloroquine diphosphate had the greatest significant inhibitory effect on autophagy in human pancreatic cancer PANC-1 cells (Fig. 2).

Figure 1

The design of the present study.

Figure 2

Autophagy inhibition caused by chloroquine diphosphate in PANC-1 cells. (A) Western blot analysis was used to detect the expression of LC3-II and P62 protein in PANC-1 cells treated with different concentrations of chloroquine diphosphate to detect the level of autophagy (*P<0.05 vs. 0 µM). (B) When the concentration of chloroquine diphosphate was 100 µM, the difference of LC3-II and P62 protein was statistically significant, compared with the control group. LC3, microtubule-associated protein-1-light-chain-3.

Quality control of samples for the ceRNA microarray and miRNA microarray chips

Electrophoresis of the RNA samples indicated that the RNA integrity number of each sample was ≥7.0 and that the 28S/18S was ≥0.7, which reached the qualifying sample standards. Thus, the samples were qualified for use in the subsequent experiments.

Assays to identify differentially-expressed mRNAs, lncRNAs, circRNAs and miRNAs

The ceRNA microarray results determined 3,966 differentially-expressed mRNAs in total. The expression levels of 2,445 mRNAs in the PANC-1 cells from the chloroquine diphosphate treatment group (autophagy suppression group) were downregulated and 1,521 mRNAs were upregulated, compared with the control group (Fig. 3). Additionally, 3,184 differentially-expressed lncRNAs were observed, including 1,637 downregulated and 1,547 upregulated (Fig. 4), whereas 9,420 circRNAs were determined to be differentially expressed, including 4,223 downregulated and 5,197 upregulated (Fig. 5). The miRNA microarray results demonstrated that the expression levels of two miRNAs (hsa-miR-663a-5p and hsa-miR-154-3p) were underexpressed in the PANC-1 cells in the autophagy-suppression group, compared with the control group. No upregulated miRNA was determined. The cutoff values of the aforementioned differentials were all set to a fold change>2 or FC<0.5 and P<0.05.

Figure 3

Scatterplots of the differentially-expressed mRNAs following autophagy being inhibited. Differentially-expressed mRNAs were screened with a microarray following autophagy being inhibited with chloroquine diphosphate at a concentration of 100 µM.

Figure 4

Scatterplots of the differentially-expressed lncRNAs following autophagy being inhibited. Differentially-expressed lncRNAs were screened with a microarray following autophagy being inhibited with chloroquine diphosphate at a concentration of 100 µM. lncRNA, long non-coding RNA.

Figure 5

Scatterplots of the differentially-expressed circRNAs following autophagy being inhibited. Differentially-expressed circRNAs were screened with a microarray following autophagy being inhibited with chloroquine diphosphate at a concentration of 100 µM. circRNA, circular RNA.

Bioinformatics analysis of the functions of the differentially-expressed mRNAs

GO and KEGG pathway analyses were performed using the DAVID database. The 3,966 differentially-expressed mRNAs in the chloroquine diphosphate treatment group (autophagy suppression group), compared with the control group, were subjected to bioinformatics analysis of their functions. The results demonstrated that these genes were concentrated in the biological processes, molecular functions and cellular components categories in the GO analysis (Table I), and were involved in the regulation of multiple KEGG signaling pathways including pathways in cancer and the mitogen-activated protein kinase signaling pathway (Table II). Of these functions, the involvement of the autophagy-associated pathway (Autophagy) was notable. Additionally, the differentially-expressed genes in this pathway included a number of recognized autophagy-associated genes, including autophagy-related 12 (ATG12), GABA type A receptor associated protein like 1 (GABARAPL1) and ULK2, which were differentially expressed in the present ceRNA microarray results also (Fig. 6). Following analyzing the clinical roles of the three genes via the Gene Expression Profiling Interactive Analysis online database, which is based on the data obtained from The Cancer Genome Atlas and GTEx, it was determined that the expression of ATG12 was upregulated in pancreatic cancer tissue, compared with non-tumor tissue. Furthermore, high expression of ATG12 was associated with significantly reduced survival time (P<0.01). Notably, the expression of ATG12, GABARAPL1 and ULK2 were increased in pancreatic cancer, compared with the control group (Fig. 6). These results validated that these autophagy-associated genes were involved in the onset and progression of pancreatic cancer.

Table I

GO analysis of differentially-expressed mRNAs following autophagy inhibition.

GO ID	Description	Type	Genes
GO:2000969	Positive regulation of α-amino-3-hydroxy-5-methyl-4-isoxazole propionate selective glutamate receptor activity	Biological_process	RELN, NLGN3, ARC and SHANK3
GO:0000727	Double-strand break repair via break-induced replication	Biological_process	CDC45, GINS2, CDC7 and GINS4
GO:0038026	Reelin-mediated signaling pathway	Biological_process	DAB2IP, RELN, LRP8 and VLDLR
GO:0003150	Muscular septum morphogenesis	Biological_process	HEY2, FZD2, VANGL2 and FZD1
GO:0050915	Sensory perception of sour taste	Biological_process	ASIC3, PKD2L1, ASIC2, PKD1L3 and ASIC1
GO:0031052	Chromosome breakage	Biological_process	HMGA2, BRCA2, BRCA1, CD19 and TP53
GO:0038171	Cannabinoid signaling pathway	Biological_process	DAGLA, CNR2, MGLL, ABHD6 and CNR1
GO:0003149	Membranous septum morphogenesis	Biological_process	ID2, WHSC1, FZD2, VANGL2 and FZD1
GO:0061156	Pulmonary artery morphogenesis	Biological_process	HEY2, STRA6, BMP4 and JAG1
GO:2001286	Regulation of caveolin-mediated endocytosis	Biological_process	UNC119, NEDD4L, CAV3 and PROM2
GO:0000942	Condensed nuclear chromosome outer kinetochore	Cellular_component	PLK1, BUB1, CCNB1 and NDC80
GO:0031262	Ndc80 complex	Cellular_component	SPC24, NUF2, SPC25 and NDC80
GO:0031298	Replication fork protection complex	Cellular_component	CDC45, MCM10, GINS2, TOP1MT and GINS4
GO:0072557	IPAF inflammasome complex	Cellular_component	CASP1, CASP12, CASP5 and NLRC4
GO:0000778	Condensed nuclear chromosome kinetochore	Cellular_component	PLK1, CENPA, NDC80, BUB1B, BUB1, CCNB1 and MIS18BP1
GO:0042555	MCM complex	Cellular_component	MCM4, TONSL, MCM5, MMS22L, MCM8, MCM2, MCMBP and MCM3
GO:0000940	Condensed chromosome outer kinetochore	Cellular_component	PLK1, SKA3, CENPF, NDC80, BUB1, BUB1B, CENPE, CCNB1 and SPDL1
GO:0005663	DNA replication factor c complex	Cellular_component	RFC4, PCNA, RFC3 and RFC5
GO:0097169	AIM2 inflammasome complex	Cellular_component	CASP1, CASP12, CASP5 and PYCARD
GO:0043240	Fanconi anemia nuclear complex	Cellular_component	FANCB, FANCC, C19ORF40, STRA13, FANCG, FANCM, FANCA and FANCE
GO:0043142	Single-stranded DNA-dependent ATPase activity	Molecular_function	PIF1, RAD18, RAD51, DNA2 and POLQ
GO:0032405	Mutlalpha complex binding	Molecular_function	TREX1, PCNA, MUTYH, MSH6 and MSH2
GO:0086006	Voltage-gated sodium channel activity involved in cardiac muscle cell action potential	Molecular_function	SCN1B, SCN4B, SCN3B and SCN2B
GO:0008574	Plus-end-directed microtubule motor activity	Molecular_function	KIF3B, KIF18B, KIF5A, KIF14, KIF5C, KIF3A, KIF20B, KIF17, KIF11, KIF26A and KIF16B
GO:0004908	Interleukin-1 receptor activity	Molecular_function	IL1R1, IL1RL2, IL1RL1, IL18R1 and IL1R2
GO:0032027	Myosin light chain binding	Molecular_function	CXCR4, IQGAP3, MYOC and MYH8
GO:0035312	5′-3′ exodeoxyribonuclease activity	Molecular_function	EXO1, DCLRE1B, EXO5 and DCLRE1A
GO:0048256	Flap endonuclease activity	Molecular_function	EXO1, MUS81, FEN1, DNA2 and FAN1
GO:0036310	Annealing helicase activity	Molecular_function	RECQL4, SMARCA1, SMARCAL1, RAD54L and BLM
GO:0097153	Cysteine-type endopeptidase activity involved in apoptotic process	Molecular_function	CASP7, CASP6, CASP12, PYCARD, CASP2, CASP8, CLC, CASP1 and CASP5

Only the top 10 signaling pathways of each classification were listed as examples. GO, Gene Ontology.

Table II

Kyoto Encyclopedia of Genes and Genomes pathway analysis of differentially-expressed mRNAs when autophagy was inhibited.

Pathway ID	Description	Genes
hsa05200	Pathways in cancer	BCR, PDGFB, WNT4, ETS1, CASP8, BRCA2, ARHGEF11, TRAF3, GNAI3 and BID
hsa04060	Cytokine-cytokine receptor interaction	CCL4, PDGFB, IL4R, CCL21, IL12A, CSF2, GH2, CCL27, ACVR1B and LEP
hsa04010	MAPK signaling pathway	PDGFB, MAP2K6, CACNG8, RPS6KA2, PLA2G4A, MAP3K4, GADD45G, CACNB3, JUN and TNFRSF1A
hsa04144	Endocytosis	FOLR3, NEDD4L, PSD4, SH3KBP1, GRK5, ERBB4, EHD2, CAV2, ITCH and EPS15L1
hsa04080	Neuroactive ligand- receptor interaction	ADRA1D, GH2, PTAFR, C5AR1, GRIN3B, HTR1B, LEP, S1PR5, LPAR2, CHRND, AGTR1, GH1, CHRM1, SSTR5, ADRA1B, GRIN3A, DRD5 and PTGER3
hsa04621	NOD-like receptor signaling pathway	CASP8, ATG12, TRAF3, TRAF5, JUN, OAS3, OAS1, MAVS, RNASEL and NLRP7
hsa00230	Purine metabolism	PFAS, PDE6G, POLD3, NME7, PDE9A, POLR3H, POLR2L, ADCY7, NME6 and NPR2
hsa04110	Cell cycle	RAD21, CDC16, E2F2, SFN, GADD45G, CCNH, CDC20, CDC25A, CCNA1 and TP53
hsa04062	Chemokine signaling pathway	CCL4, CCL21, GNAI3, PIK3R5, GRK5, CCL27, PIK3CB, SHC1, VAV1 and CCL4L1
hsa04210	Apoptosis	CASP8, BID, PIK3R5, CTSC, LMNB1, PIK3CB, EIF2AK3, GADD45G, CASP6 and TRADD

Only 10 genes were listed as examples in each pathway. MAPK, mitogen-activated protein kinase.

Figure 6

Expression level, stage and overall survival curves of ATG12, GABARAPL1 and ULK2 in patients with pancreatic cancer. (A) The fold change of ATG12 (fold change=2.34), GABARAPL1 (fold change=3.90) and ULK2 (fold change=2.29) in the chloroquine diphosphate treatment group, compared with the control group. (B) The expression of ATG12 was upregulated in pancreatic cancer tissues, compared with normal tissues [T, n=179; N, n=171]. (C) ATG12 expression in different stages of pancreatic cancer is not significantly different (P=0.6140). (D) Patients with pancreatic cancer and a low expression of ATG12 have an increased survival time, compared with those with a high expression of ATG12 (HR=1.8; P=0.0042). (E) GABARAPL1 has an increased expression in pancreatic cancer tissues, compared with normal tissues [T, n=179; N, n=171]. (F) The association between GABARAPL1 and stage of patients with pancreatic cancer (P=0.5940). (G) The effect of GABARAPL1 expression on the overall survival time of patients with pancreatic cancer is not significant (HR=0.72; P=0.1100). (H) ULK2 has an increased expression in pancreatic cancer tissues, compared with normal tissues [T, n=179; N, n=171]. (I) The association between ULK2 and stage of patients with pancreatic cancer (P=0.0248). (J) The effect of ULK2 expression on the overall survival time of patients with pancreatic cancer is not significant (HR=0.79; P=0.2500). ATG12, autophagy-related 12; GABARAPL1, GABA type A receptor associated protein like 1; ULK2, Unc-51 like autophagy activating kinase 2; HR, hazard ratio; TPM, transcript per million.

Prediction of target genes of the differentially-expressed miRNAs and the bioinformatics analysis of their functions

These results were integrated with the differentially-expressed mRNAs that were actually measured. The co-expressed genes were subjected to GO and KEGG pathway analyses using the DAVID database.

Prediction of miR-663a-5p target genes and bioinformatics analysis of their functions

A total of 1,726 target genes of miR-663a-5p were detected by at least 6 platforms using the aforementioned 12 online miRNA target gene prediction databases. Following cross-checking with the differentially-expressed genes in the autophagy suppression group derived from chloroquine diphosphate treatment, 462 co-expressed genes were obtained (Fig. 7). These genes were subjected to GO and KEGG pathway analyses. The results demonstrated that these genes were concentrated in the biological processes, molecular functions and cellular components categories in the GO analysis (Fig. 8A, Table III) and were involved in the regulation of multiple KEGG signaling pathways (Table III), including the aldosterone-regulated sodium reabsorption and Wnt signaling pathways (Fig. 8B-D). However, there were no significant pathways determined by PANTHER analysis (P>0.05).

Figure 7

Potential target genes of microRNA-663a-5p associated with autophagy. The 462 genes from target prediction and differentially-expressed mRNAs were depicted in Compendia expression profiles.

Figure 8

Circular visualization of the results of GO and KEGG pathway analysis of miR-663a-5p target genes. (A) The horizontal axis depicts the number of downregulated genes of miR-663a-5p. The vertical axis depicts the GO categories. (B) Circle plot demonstrating the significance of the pathway enriched by the target genes of miR-663a-5p. Inner plot color corresponds to the z-score. The outer ring displays scatterplots of the expression levels (logFC) for the genes in each term. (C) Cluster plot depicting a circular dendrogram of the target genes of miR-663a-5p. The inner ring depicts the color-coded logFC, the outer ring depicts the assigned functional terms. (D) Chord plot depicting the target genes of miR-663a-5p associations via ribbons to their assigned pathways. Colored rectangles represent the logFC of genes. GO, Gene Ontology; miR, microRNA; logFC, log fold change; PI3K, phosphoinositide 3-kinase.

Table III

GO and KEGG analysis of potential target genes of microRNA-663a-5p related to autophagy.

Term	Description	Type	Count	P-value	Genes
GO:0048562	Embryonic organ morphogenesis	Biological process	11	1.66×10−3	DLX2, FOXL2, TBX15, MAFB, CHST11, NKX3-2, SOBP, TBX1, ZEB1 and PAX2
GO:0045449	Regulation of transcription	Biological process	87	1.71×10−3	ZNF451, NR6A1, ZNF250, CBX2, ZEB1, MED22, PAX2, CBFA2T3, ZNF345 and TGFB1
GO:0006350	Transcription	Biological process	73	1.78×10−3	NR6A1, ZNF451, ZNF250, CBX2, MED22, ZEB1, PAX2, CBFA2T3, ZNF345 and CRX
GO:0043583	Ear development	Biological process	9	2.35×10−3	KCNMA1, MAFB, NKX3-2, SOBP, TBX1, PAX2, SLC9A3R2, TGFB1 and CDH23
GO:0017015	Regulation of transforming growth factor β receptor signaling pathway	Biological process	6	2.54×10−3	DAND5, CHST11, TGFB1I1, ZEB1, PRDM16 and TGFB1
GO:0048839	Inner ear development	Biological process	8	3.27×10−3	KCNMA1, MAFB, SOBP, TBX1, PAX2, SLC9A3R2, TGFB1 and CDH23
GO:0048568	Embryonic organ development	Biological process	12	3.56×10−3	DLX2, FOXL2, TBX15, MAFB, FOXF1, CHST11, NKX3-2, SOBP, TBX1 and ZEB1
GO:0006357	Regulation of transcription from RNA polymerase II promoter	Biological process	31	3.82×10−3	TADA3, CRTC1, NR6A1, CBX2, ZEB1, MED22, PAX2, PRDM16, ZNF345 and TGFB1
GO:0006355	Regulation of transcription DNA-dependent	Biological process	62	3.88×10−3	NR6A1, ZNF250, CBX2, MED22, ZEB1, CBFA2T3, PAX2, ZNF345, TGFB1 and CRX
GO:0048598	Embryonic morphogenesis	Biological process	17	3.96×10−3	FOXL2, TBX15, MAFB, ARFRP1, TP53, SOBP, TBX1, ZEB1, CELSR1 and PAX2
GO:0030054	Cell junction	Cellular component	27	1.02×10−4	ACHE, CLDN9, PANX2, CLDN6, GRIK5, ZNRF1, ITSN1, SYNGR1, CALB2 and RIMS3
GO:0045202	Synapse	Cellular component	21	1.45×10−4	KCNMA1, DLGAP1, ACHE, CDK5R1, ARC, EFNA2, SYT11, GRIK5, BCAN and BSN
GO:0044459	Plasma membrane part	Cellular component	68	4.22×10−3	KCNC1, PCDHA8, CLDN9, CD8A, CLDN6, GRIK5, ANPEP, SYNGR1, ITSN1 and ZNRF1
GO:0019717	Synaptosome	Cellular component	7	1.19×10−2	PVRL1, SNPH, DLG4, BSN, LGI3, ITSN1 and RNF40
GO:0030173	Integral to golgi membrane	Cellular component	5	1.53×10−2	ST6GAL1, ST8SIA4, RER1, STEAP2 and ABO
GO:0031228	Intrinsic to golgi membrane	Cellular component	5	1.92×10−2	ST6GAL1, ST8SIA4, RER1, STEAP2 and ABO
GO:0005794	Golgi apparatus	Cellular component	29	3.41×10−2	SLC9A8, ACHE, APC2, C6ORF25, PPIL2, ARFRP1, RER1, VPS53, CBFA2T3 and NUFIP2
GO:0005886	Plasma membrane	Cellular component	99	4.57×10−2	KCNC1, ALPPL2, PLXNA1, ATP1B2, C6ORF25, EFNA2, FGFRL1, GRIK5, ZNRF1, TLR6 and ITSN1
GO:0031012	Extracellular matrix	Cellular component	14	4.66×10−2	ACHE, ADAMTS14, ADAMTSL2, ADAMTS15, OLFML2A, BCAN, COL5A1, MMP25, TGFB1 and WNT7B
GO:0005667	Transcription factor complex	Cellular component	10	4.72×10−2	BRF1, TADA3, TBX2, MAFB, FOXF1, NR6A1, TP53, NPAS4, ZEB1 and CRX
GO:0046872	Metal ion binding	Molecular function	123	1.21×10−2	KCNC1, SLC9A8, ALAD, FSTL4, ALPPL2, ATP1B2, ZNF451, ZNF250, RNF216 and ITSN1
GO:0008270	Zinc ion binding	Molecular function	74	1.44×10−2	ALPPL2, ALAD, ZCCHC24, ZNF451, NR6A1, SOBP, ZNF250, ANPEP, RNF216 and ZEB1
GO:0043169	Cation binding	Molecular function	123	1.64×10−2	KCNC1, SLC9A8, ALAD, ALPPL2, ATP1B2, ZNF451, FSTL4, ZNF250, RNF216 and ITSN1
GO:0019209	Kinase activator activity	Molecular function	4	2.32×10−2	CDK5R1, MADD, ITSN1 and TAB1
GO:0043167	Ion binding	Molecular function	123	2.58×10−2	KCNC1, SLC9A8, ALAD, ALPPL2, ATP1B2, ZNF451, FSTL4, ZNF250, RNF216 and ITSN1
GO:0003677	DNA binding	Molecular function	72	3.34×10−2	PRR12, NR6A1, ZNF451, ZNF250, CBX2, ZEB1, PAX2, CBFA2T3, ZNF345 and CRX
GO:0046914	Transition metal ion binding	Molecular function	84	3.45×10−2	STEAP3, ALPPL2, ALAD, ZCCHC24, ZNF451, NR6A1, SOBP, ZNF250, ANPEP and RNF216
GO:0030955	Potassium ion binding	Molecular function	8	4.02×10−2	KCNMA1, KCNC1, SLC12A7, PDXK, ATP1B2, KCNK5, HCN4 and KCNG1
GO:0030528	Transcription regulator activity	Molecular function	49	4.27×10−2	NR6A1, CBX2, ZEB1, MED22, CBFA2T3, ZNF345, TGFB1, CRX, FOXF1 and ZNF445
GO:0016419	S-malonyltransferase activity	Molecular function	2	4.90×10−2	MCAT and FASN
hsa04960	Aldosterone-regulated sodium reabsorption	KEGG pathway	5	9.51×10−3	ATP1B2, PIK3R5, NEDD4L, SFN and SLC9A3R2
hsa04310	Wnt signaling pathway	KEGG pathway	8	2.92×10−2	WNT7B, SOST, VANGL1, APC2, TP53, FZD1, FRAT1 and FZD4
hsa05217	Basal cell carcinoma	KEGG pathway	5	3.03×10−2	WNT7B, APC2, TP53, FZD1 and FZD4
hsa00061	Fatty acid biosynthesis	KEGG pathway	3	3.08×10−2	MCAT, FASN and ACSL6
hsa04151	PI3K-Akt signaling pathway	KEGG pathway	14	3.48×10−2	FGF18, IL2RB, LAMB3, FLT1, IL4R, CSF1, EFNA2, TP53, TNN and EFNA5
hsa05205	Proteoglycans in cancer	KEGG pathway	9	6.69E×10−2	WNT7B, CTTN, EZR, ANK1, TP53, FZD1, PIK3R5, FZD4 and TGFB1
hsa04724	Glutamatergic synapse	KEGG pathway	6	9.94×10−2	DLGAP1, SLC1A7, DLG4, GRIK5, CACNA1D and SHANK3

Only 10 genes were listed as examples in each pathway. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PI3K, phosphoinositide 3-kinase.

Prediction of miR-154-3p target genes and bioinformatics analysis of their functions

A total of 294 target genes of miR-154-3p were detected by at least 6 programs using the aforementioned 12 online miRNA target gene prediction databases. Following cross-checking with the differentially-expressed genes in the autophagy suppression group derived from chloroquine diphosphate treatment, 294 co-expressed genes were obtained (Fig. 9). These genes were subjected to GO and KEGG pathway analyses. The results demonstrated that these genes were concentrated in the biological processes, molecular functions and cellular components categories in the GO analysis (Fig. 10A, Table IV), and were involved in the insulin signaling, Forkhead Box O (FoxO) signaling and proteoglycans in cancer pathways (Table IV, Fig. 10B-D). Additionally, the results of PANTHER analysis revealed these genes were concentrated in a numbered of pathways, including the p53 pathway by glucose deprivation and the classical gastrin cholecystokinin receptor (CCKR) signaling map (Table V, Fig. 11).

Figure 9

Potential target genes of microRNA-154-3p associated with autophagy. The 294 genes from target prediction and differentially-expressed mRNAs were depicted in Compendia expression profiles.

Figure 10

Circular visualization of the results of GO and KEGG pathway analysis of miR-154-3p target genes. (A) The horizontal axis depicts the number of downregulated genes of miR-154-3p. The vertical axis depicts the GO categories. (B) Circle plot depicting the significance of the pathway enriched by the target genes of mir-154-3p. Inner plot color corresponds to the z-score. The outer ring displaying scatterplots of the expression levels (logFC) for the genes in each term. (C) Cluster plot displaying a circular dendrogram of the target genes of mir-154-3p. The inner ring depicting the color-coded logFC, the outer ring depicting the assigned functional terms. (D) Chord plot depicting the target genes of miR-154-3p associations via ribbons to their assigned pathways. Colored rectangles represent the logFC of genes. miR, microRNA; logFC, log fold change; FoxO, Forkhead box O; MAPK, mitogen-activated protein kinase.

Table IV

GO and KEGG analysis of potential target genes of microRNA-154-3p associated with autophagy.

Term	Description	Type	Count	P value	Genes
GO:0007605	Sensory perception of sound	Biological process	11	1.67×10−3	NAV2, ASIC2, DCDC2, SOBP, TBX1, TIMM13, CACNA1D, FZD4, SLC52A3 and DNM1
GO:0030512	Negative regulation of transforming growth factor β receptor signaling pathway	Biological process	7	4.81×10−3	DAND5, ADAMTSL2, ZNF451, CHST11, TGFB1I1, PRDM16 and TGFB1
GO:0072207	Metanephric epithelium development	Biological process	3	5.74×10−3	WNT7B, OSR1 and PAX2
GO:0070507	Regulation of microtubule cytoskeleton organization	Biological process	4	6.52×10−3	DIXDC1, CDK5R1, ATAT1 and EFNA5
GO:0045216	Cell-cell junction organization	Biological process	4	9.16×10−3	CLDN9, CLDN6, MARVELD3 and TGFB1
GO:0008285	Negative regulation of cell proliferation	Biological process	19	9.45×10−3	TP53I11, NACC2, CLMN, KLF10, FGFRL1, TP53, GPER1, ZEB1, CBFA2T3 and TGFB1
GO:0031497	Chromatin assembly	Biological process	3	1.53×10−2	CDAN1, TP53 and CHAF1A
GO:0006366	Transcription from RNA polymerase II promoter	Biological process	22	1.62×10−2	MAF, MAFF, FOXL2, POLR2L, MAFB, SOX12, TP53, ARID3B, FOSB and NPAS4
GO:0008104	Protein localization	Biological process	6	1.69×10−2	TNS4, MALL, CDAN1, TP53, GRASP and AKAP3
GO:0000122	Negative regulation of transcription from RNA polymerase II promoter	Biological process	28	1.98×10−2	NR6A1, CBX2, ZEB1, PRDM16, ZNF345, TGFB1, KANK2, AHRR, EZR and OSR1
GO:0030054	Cell junction	Cellular component	24	6.45×10−4	SH3PXD2B, ARC, DLGAP1, ACHE, SH3PXD2A, CACNG8, SYT11, GRIK5, LRRC4B and BSN
GO:0014069	Postsynaptic density	Cellular component	13	1.51×10−3	DLGAP1, ARC, CDK5R1, CACNG8, SYT11, GRIK5, BSN, GPER1, SHANK3 and CNIH2
GO:0005737	Cytoplasm	Cellular component	152	2.33×10−3	RBPMS2, STIL, PLXNA1, ATP1B2, LRRC4B, MED22, CALB2, TGFB1, KANK2 and PACSIN1
GO:0043197	Dendritic spine	Cellular component	9	2.60×10−3	CDK5R1, ARC, CTTN, CNIH2, SYT11, DLG4, ASIC2, SEZ6 and SHANK3
GO:0044297	Cell body	Cellular component	7	3.72×10−3	GNAZ, EZR, SYT11, PACRG, FADD, DISC1 and EPO
GO:0005634	Nucleus	Cellular component	151	1.30×10−2	ALAD, PLXNA1, ZNF451, ZNF250, CBX2, RNF216, ITPKB, CALB2, TGFB1 and ANK1
GO:0030424	Axon	Cellular component	12	1.72×10−2	CDK5R1, ATAT1, STMN3, SYT11, FKBP15, LMTK3, BSN, GPER1, LDLRAP1 and KIF21B
GO:0005884	Actin filament	Cellular component	6	1.89×10−2	CTTN, EZR, APC2, AIF1L, FKBP15 and MYO9B
GO:0005886	Plasma membrane	Cellular component	115	3.45×10−2	KCNC1, SLC9A8, ALPPL2, PLXNA1, ATP1B2, EFNA2, C6ORF25, FGFRL1, GRIK5 and TLR6
GO:0048471	Perinuclear region of cytoplasm	Cellular component	23	4.02×10−2	CDK5R1, ACHE, STC2, APC2, CSF1, BRSK2, VPS53, MYO9B, GPER1 and NDOR1
GO:0001228	Transcriptional activator, activity RNA polymerase II transcription regulatory region sequence-specific binding	Molecular function	8	8.50×10−3	MAF, MAFF, DLX2, MAFB, FOXF1, TP53, NEUROD2 and ARID3B
GO:0046872	Metal ion binding	Molecular function	68	8.68×10−3	GNAZ, STEAP3, SGSH, ALPPL2, ALAD, ZNF451, SOBP, ZNF250, RNF216 and ZEB1
GO:0000978	RNA polymerase II core promoter proximal region sequence-specific DNA binding	Molecular function	17	1.37×10−2	FOXL2, TBX15, NACC2, TBX2, MAFB, NR6A1, ARID3B, FOSB, NPAS4 and DDN
GO:0048495	Roundabout binding	Molecular function	3	1.90×10−2	MYO9B, TGFB1I1 and TPBGL
GO:0003779	Actin binding	Molecular function	14	1.90×10−2	KCNMA1, DIXDC1, FMNL3, CLMN, MYO9B, SHANK3, TNS4, EPB41L1, EZR and PACRG
GO:0030295	Protein kinase activator activity	Molecular function	4	2.46×10−2	DUSP19, MADD, RPTOR and EPO
GO:0098641	Cadherin binding involved in cell-cell adhesion	Molecular function	14	2.58×10−2	EPS15L1, SFN, SLC9A3R2, TNKS1BP1, PAK6, CTTN, EPB41L1, EZR, FASN and SPTBN1
GO:0016874	Ligase activity	Molecular function	13	3.31×10−2	DTX3L, PPIL2, ZNF451, RNF216, ZNRF1, TTLL11, RNF165, NEURL1B, PELI2 and NEDD4L
GO:0030276	Clathrin binding	Molecular function	5	3.75×10−2	TOM1L2, SYT11, RPH3AL, C2CD4C and LDLRAP1
GO:0004314	[acyl-carrier-protein] S-malonyltransferase activity	Molecular function	2	4.81×10−2	MCAT and FASN
hsa04910	Insulin signaling pathway	KEGG pathway	8	2.79×10−3	CBLB, SOS1, PRKAB2, MAPK8, RPS6KB1, SOCS4, PRKACB and PPARGC1A
hsa04068	FoxO signaling pathway	KEGG pathway	7	1.01×10−2	NLK, SOS1, PRKAB2, MDM2, FBXO32, MAPK8 and FOXO3
hsa04931	Insulin resistance	KEGG pathway	6	1.62×10−2	RPS6KA2, PRKAB2, MAPK8, RPS6KB1, OGT and PPARGC1A
hsa05205	Proteoglycans in cancer	KEGG pathway	8	1.98×10−2	CBLB, CTTN, ERBB4, SOS1, MDM2, RPS6KB1, PRKACB and FZD4
hsa04012	ErbB signaling pathway	KEGG pathway	5	3.14×10−2	CBLB, ERBB4, SOS1, MAPK8 and RPS6KB1
hsa04010	MAPK signaling pathway	KEGG pathway	8	6.08×10−2	RPS6KA2, ARRB1, NLK, SOS1, FGF13, MAPK8, PRKACB and MAP3K13
hsa00410	β-Alanine metabolism	KEGG pathway	3	6.72×10−2	ALDH6A1, ALDH1B1 and EHHADH
hsa00512	Mucin type O-Glycan biosynthesis	KEGG pathway	3	6.72×10−2	GALNT4, GALNT18 and B4GALT5
hsa04920	Adipocytokine signaling pathway	KEGG pathway	4	7.13×10−2	PRKAB2, MAPK8, PPARGC1A and CAMKK2
hsa04917	Prolactin signaling pathway	KEGG pathway	4	7.13×10−2	SOS1, MAPK8, SOCS4 and FOXO3

Only 10 genes were listed as examples in each pathway. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MAPK, mitogen-activated protein kinase; FoxO, Forkhead box O.

Table V

PANTHER pathway analysis of potential target genes of miR-663a-5p and miR-154-3p in pancreatic cancer.

miRNA	Geneset	Description	Count	P-value	E	Overlapping genes	URL
miR-663a-5p	P00037	Ionotropic glutamate receptor pathway	3	0.0852	1.05	GRIK5, SLC1A7 and SHANK3	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P00037
	P02743	Formyltetrahydroformate biosynthesis	1	0.1319	0.14	MTR	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P02743
	P04392	P53 pathway feedback loops 1	1	0.1319	0.14	TP53	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P04392
	P00012	Cadherin signaling pathway	6	0.1412	3.56	PCDHA8, CDH23, WNT7B, FZD1, FZD4 and CELSR1	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P00012
	P05726	2-arachidonoylglycerol biosynthesis	1	0.1522	0.16	DAGLA	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P05726
	P00020	FAS signaling pathway	2	0.1614	0.72	DFFB and FADD	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P00020
	P00029	Huntington disease	5	0.1702	2.96	ARL4C, DLG4, GRIK5, PACSIN1 and TP53	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P00029
	P00057	Wnt signaling pathway	9	0.2380	6.84	APC2, PCDHA8, CDH23, TP53, WNT7B, FZD1, FZD4, KREMEN1 and CELSR1	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P00057
	P00045	Notch signaling pathway	2	0.2383	0.93	HEYL and HES2	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P00045
	P02746	Heme biosynthesis	1	0.2467	0.28	ALAD	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P02746
miR-154-3p	P00032	Insulin/IGF pathway-MAPK kinase/MAPK cascade	3	0.0116	0.48	RPS6KA2, RPS6KB1 and SOS1	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P00032
	P06959	CCKR signaling map	7	0.0231	2.92	FOXO3, CXCL2, PRKACB, MAPK8, RPS6KB1, SOS1 and TAC1	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P06959
	P04372	5-Hydroxytryptamine degradation	2	0.0351	0.30	ALDH1L2 and ALDH1B1	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P04372
	P00018	EGF receptor signaling pathway	5	0.0398	1.92	ERBB4, MAPK8, SOS1, CBLB and PHLDB2	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P00018
	P00048	PI3K pathway	3	0.0420	0.78	FOXO3, RPS6KB1 and SOS1	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P00048
	P04397	p53 pathway by glucose deprivation	2	0.0467	0.35	PRKAB2 and RPS6KB1	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P04397
	P00041	Metabotropic glutamate receptor group I pathway	2	0.0594	0.40	GRIN2A and PRKACB	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P00041
	P00001	Adrenaline and noradrenaline biosynthesis	2	0.0733	0.45	SLC6A19 and SLC6A15	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P00001
	P00035	Interferon-γ signaling pathway	2	0.0781	0.47	SOCS4 and MAPK8	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P00035
	P00058	mRNA splicing	1	0.0808	0.08	SNRNP40	http://www.pantherdb.org/pathway/pathwayDiagram.jsp?catAccession=P00058

miR, microRNA; MAPK, mitogen-activated protein kinase; PI3K, phosphoinositide 3-kinase; IGF, interferon growth factor; CCKR, classical gastrin cholecystokinin receptor.

Figure 11

Bubble Chart of the results of PANTHER analyses. (A) Pathway visualization of the targeted genes of miR-154-3p. (B) Pathway visualization of the targeted genes of miR-154-3p. miR, microRNA; CCKR, classical gastrin cholecystokinin receptor; MAPK, mitogen-activated protein kinase; MAP2K, mitogen-activated protein kinase kinase; EGF, epidermal growth factor; PI3K, phosphoinositide 3-kinase.

Prediction of ceRNAs of the differentially-expressed miRNAs

Different online platforms were used to predict the target genes of the differentially-expressed miRNAs. The target circRNAs and lncRNAs of the differentially-expressed miRNAs were searched based on the prediction of the target miRNAs of the differentially-expressed circRNAs and lncRNAs revealed by microarray detection, and the differentially-expressed miRNAs that were detected. The ceRNAs that were negatively associated with the expression of the miRNAs were selected.

Prediction of ceRNAs of miR-663a-5p

A total of 21 differentially-expressed circRNAs that targeted miR-663a-5p were determined, of which nine (hsa_circ_0003176, hsa_circ_0048579, hsa_circ_0063706, hsa_circ_0071922, hsa_circ_0078989, hsa_circ_0079319, hsa_circ_0083080, hsa_circ_0089643 and hsa_circ_0090372) indicated negative associations with miR-663a-5p expression. hsa_circ_0071922 had the highest number of binding sites for miR-663a-5p, with six. Simultaneously, 45 differentially-expressed lncRNAs that targeted miR-663a-5p were determined, of which 8 lncRNAs (RP11-59C5.3, RP13-516M14.8, RP11-196G18.24, AJ006995.3, AC024560.2, PPP1R1C, LINC00595 and HAGLROS) were negatively associated with miR-663a-5p expression. All of these lncRNAs had one miR-663a-5p binding site (Table VI). Subsequently, 144 common predicted targets from at least 8 among 12 programs were selected to overlap with the 3,966 differentially-expressed genes following autophagy inhibition. Collectively, 46 potential targets were determined. Thus, a ceRNA hypothesis figure with miR-663a-5p, 9 circRNAs, 8 lncRNAs and 46 genes was depicted in Fig. 12.

Table VI

ceRNAs and their sequences of miR-663a-5p and miR-154-3p.

miRNA	ceRNA	Sequences
miR-663a-5p	hsa_circ_0003176	5′-ACATGGAGCTGCACAGAATGTCAAGAATGCAAAGAGCAGAACCTTGCGCGGGCACAGGGC-3′
	hsa_circ_0048579	5′-TTCCTTAGGTTACAAAACAAAACAGGGAGAGAAAAAGCAAGTTCTCCGGGGCTTGGCGTA-3′
	hsa_circ_0063706	5′-ACTGGTCAAGTGAGTGACCTGCAGGAAGAAGGCAAGAATGCCATCAACTCACCGATGTCC-3′
	hsa_circ_0071922	5′-GACCTCTGTTGTACCTGTAATAAATATATAGAAAAAGCAAACACATACACCAATGCGGAT-3′
	hsa_circ_0078989	5′-TGTGGATGCTGTCAAACTGGGTCTACCGATTCTTCGCTGCTGCTGCCTTACCGCCGAGAA-3′
	hsa_circ_0079319	5′-TTGTCAAGTGTCAATAAAAGCATCATGTAATTTACCGAACGCCCGCAGCAGGGTCAGAAG-3′
	hsa_circ_0083080	5′-TGCATTCATTTTTTACTTTTATTAAAGGTTCAAAACCAACACGCCCCTGGGCGGTGCGCG-3′
	hsa_circ_0089643	5′-CGTCACGGCGCGCTACTACGGGAAGTCGCTGCCGTTCGGTGCGCAGTCCACGCAGCGCGG-3′
	hsa_circ_0090372	5′-AGCACCTAATTAAAGGGGCTGAAAGTCTGAGCTCTGGCCGGCCCCGGCGATTGGTCACCG-3′
	RP11-59C5.3	5′-CTTGTTTGTGCCAGGGAAAGATGGAAAATTTTAGTTTATATTAAACACGCCTTTGAAGAC-3′
	RP13-516M14.8	5′-AAACACCTGTCAGATTTCTAAGACCTGAAGCCTTGAGAAAGCATCGAAAAGTATTCAAGT-3′
	RP11-196G18.24	5′-CGCTACATAGTACGTATTATGACCAATGGATATCATTCAGAGTTCTAGTTTACATATGCT-3′
	AJ006995.3	5′-GGGAACAATTTATCTATGTTAAAAAGATGTGTGTGTGTGAATGGATGTATATACGTGTGC-3′
	AC024560.2	5′-TCACCGTGAAATATCTTCTGAGACACAACAATCTCCCTGAATTAAAACTACAGACTCTAA-3′
	PPP1R1C	5′-ACTTTTCTGAGTATACCATGGAATTCCACTGCTTGACTTCCAGAAGCATCCTCCATCTCT-3′
	LINC00595	5′-GCTGGCATCTAACGGCTTTTGTTAATTACTATATCCCAAAGAAAATGAAATAAAGGAGGG-3′
	HAGLROS	5′-GAAGTAGGCTAAGACTGCTGTTGAACTTGAAGCTAAACTTGAATGCCATTTAAAATAACA-3′
miR-154-3p	hsa_circ_0000156	5′-TGCTTTCAAGAATCATAGTGAAAGATTCGTTTTCATTGCAGAGTGGTATGATCCAAATGC-3′
	hsa_circ_0004089	5′-CACCTGCAGTAGGAAAAAGAAGCAACTAACAAAACACTGTGATAATAAGGATTATTCAGT-3′
	hsa_circ_0006461	5′-GTATGATGACATCATCCCCCCGAAGGAGTGGAAGCCGCGGCAGACGTATGATGACATCGA-3′
	hsa_circ_0015157	5′-GGAAAGAAAAAATCATAGTGAAAGATTCGTTTTCATTGCAGAGTGGTATGATCCAAATGC-3′
	hsa_circ_0038665	5′-TAAAAAATGTTTGGAGTGTAAATATACTGATGATAATGCTCTGATTCCTAAGAATTCTTC-3′
	RP11-686O6.1	5′-CCCTGCTTATGATGTAAAAATGGAGGCAGAAGACCTAAATCTAAGGCAATATACAATTTA-3′
	LINC01140	5′-CAGATGGGGTTGAAAACTCCCCTCAAAATGCAGTCTTGTTAAATTAAAAACTTGATCCAG-3′

miR, microRNA; ceRNA, competitive endogenous RNA.

Figure 12

ceRNA network of miR-663a-5p. The predicted ceRNA of miR-663a-5p included 9 upregulated circRNAs, 8 upregulated lncRNAs, miR-663a-5p and 46 differentially-expressed mRNA, which were differentially expressed in the ceRNA microarray results as well. Nodes in green, yellow, red and blue represent circRNAs, lncRNAs, miRNA-663a-5p and mRNAs, respectively. miR, microRNA; circRNA, circular RNA; lncRNA, long non-coding RNA; ceRNA, competitive endogenous RNA.

Prediction of ceRNAs of miR-154-3p

A total of 6 differentially-expressed circRNAs that targeted miR-154-3p were determined, of which five (hsa_circ_0000156, hsa_circ_0004089, hsa_circ_0006461, hsa_circ_0015157 and hsa_circ_0038665) indicated a negative association with miR-154-3p expression. All of these circRNAs had two miR-154-3p binding sites. Simultaneously, 16 differentially-expressed lncRNAs that targeted miR-154-3p were determined, of which two (RP11-686O6.1 and LINC01140) were negatively associated with miR-154-3p expression. These lncRNAs had one miR-154-3p binding site (Table VI). Similarly to miR-663a-5p, 67 common predicted targets of miR-154-3p from at least 8 among 12 platforms were collected, and then they were intersected into the 3,966 differentially-expressed genes following autophagy suppression. Eventually, 11 potential targets were collected. Hence, a ceRNA hypothesis network with miR-154-3p, 5 circRNAs, 2 lncRNAs and 11 genes was presented (Fig. 13).

Figure 13

ceRNA network of miR-154-3p. The predicted ceRNA of miR-154-3p included 5 upregulated circRNAs, 2 upregulated lncRNAs, miR-154-3p and 11 differentially-expressed mRNA, which were differentially expressed in the ceRNA microarray results as well. Nodes in yellow, green, red and blue represent circRNAs, lncRNAs, miRNA-154-3p and mRNAs, respectively. The arrows indicate the direction of transcription. miR, microRNA; circRNA, circular RNA; lncRNA, long non-coding RNA; ceRNA, competitive endogenous RNA.

Meta-analysis of miR-663a-5p and miR-154-3p expression based on RT-qPCR and GEO data

A comprehensive analysis for the expression of miR-663a-5p and miR-154-3p in pancreatic cancer tissues was performed based on RT-qPCR and GEO data. The combined standardized mean difference (SMD) values of miR-663a-5p was −0.203 [95% confidence interval (CI), −0.675-0.269; Fig. 14A)], which indicated that miR-663a had low expression in pancreatic cancer tissues; however, heterogeneity existed (I2=63.6%; P=0.011; Fig. 14B). Sensitivity analysis and publication bias of miR-663a-5p were depicted in Fig. 14C. The combined SMD values of miR-154-3p was −0.434 (95% CI, −1.079-0.212; Fig. 15A)], which indicated that miR-154-3p had low expression in pancreatic cancer tissues; however, heterogeneity existed (I2=76.0%; P=0.001; Fig. 15B). Sensitivity analysis and publication bias of miR-154-3p were depicted in Fig. 15C.

Figure 14

Meta-analysis of miR-663a-5p in pancreatic cancer. (A) A forest plot of miR-663a-5p expression. (B) Sensitivity analysis for miR-663a-5p and (C) Begg’s funnel plot of miR-663a-5p expression (Pr>|z|=0.548). miR, microRNA; SMD, standardized mean difference; s.e., standard error of the effect size.

Figure 15

Meta-analysis of miR-154-3p in pancreatic cancer. (A) A forest plot of miR-154-3p expression. (B) Sensitivity analysis for miR-154-3p. (C) Begg’s funnel plot of miR-154-3p expression (Pr>|z|=0.707). miR, microRNA; SMD, standardized mean difference; s.e., standard error of the effect size.

Discussion

Pancreatic cancer is one of the common malignancy types of the digestive system. Due to the lack of an effective early diagnosis, numerous patients are already in the advanced stage of the cancer when diagnosed (85). Currently, pancreatic cancer treatment remains dominated by surgical resection, which has a low five-year survival rate (86-88). Therefore, there is an urgent requirement for a series of effective markers for pancreatic cancer to change the current status of the poor efficacy of individualized treatment. Recent studies demonstrated that although ncRNAs do not encode proteins, they serve a pivotal role in the regulation of a variety of malignant tumor types, including gastric (89), pancreatic (90), prostate (91) and breast cancer (92). Over the past two years, the extensively investigated ncRNAs, including lncRNAs and miRNAs, and circRNAs, which were newly revealed to have similar regulatory functions, particularly to miRNAs, were indicated to also serve crucial roles (44,93-97). Abnormal expression of a large number of miRNAs has been determined in malignant tumor tissues (98-100). For example, it has been reported that miR-221/222 and miR-15b serve a role in causing malignant tumors (101). Additionally, it was determined that miR-451 and miR-126 have an abnormal expression in lung cancer (102), and the targeted delivery of miRNA therapeutics was a promising strategy for cancer (103). This abnormal expression pattern affects the occurrence, development and prognosis of the tumor by directly regulating the biological functions of the targeted mRNAs (102,104-112). Studies have demonstrated that a variety of miRNAs exhibit abnormal expression patterns in pancreatic cancer cells that affect the occurrence, development and prognosis of pancreatic cancer; for example, the upregulated expression of miR-10b, miR-21, miR23a and miR-27a in pancreatic cancer affects cell growth, proliferation and apoptotic metastasis by targeting programmed cell death 4, BTG anti-proliferation factor 2, neural precursor cell expressed, developmentally downregulated 4-like, phosphatase and tensin homolog, HIV-1 Tat interactive protein 2 and p16 (113-117).

With the increased understanding of the regulatory mechanism of ncRNAs, ceRNAs have attracted increasing attention as important regulators of miRNA activity. ceRNA refers to lncRNAs and circRNAs that suppress miRNA expression through targeted binding to the miRNA, and thus regulate its activity, which affects the development, progression and prognosis of tumors (118-120). Studies have indicated that lncRNA-metastasis associated lung adenocarcinoma transcript 1 can target miR-124 to reduce its expression and activate cyclin-dependent kinase 4, thereby accelerating the progression of breast cancer (121-123). A recent study on pancreatic cancer demonstrated that lncRNA-urothelial cancer associated 1 was overexpressed in pancreatic cancer tissues and reduced miR-135a expression by adsorbing it, thereby inhibiting genes associated with tumor growth and metastasis (124). These examples highlight the significance of the interactions among lncRNAs, miRNAs and mRNAs for cancer diagnosis and treatment.

Therefore, our aim was to mine pancreatic cancer autophagy-associated ceRNAs and miRNAs through assays with ceRNA and miRNA microarrays, respectively, to reveal their targets using bioinformatics methods, construct ceRNA, miRNA and mRNA pathways, predict the potential molecular mechanisms of the ceRNA, miRNA and mRNA pathways in the autophagy of pancreatic cancer, and provide novel ideas and directions for further studies of pancreatic cancer.

Firstly, the differentially-expressed mRNAs were analyzed. The results demonstrated that PANC-1 cells treated with chloroquine diphosphate had multiple differentially-expressed mRNAs, compared with the control group, using the FC value FC>2 or FC<0.5 as the threshold. The bioinformatics analysis of these genes indicated that as the autophagic level changed, the differentially-expressed genes were primarily concentrated in tumor-associated and pancreatic cancer-associated pathways. Notably, the autophagy-associated pathway, in which three genes (ATG12, GABARAPL1 and ULK2) were abnormally expressed when the autophagic level, was inhibited. ATG12 is an autophagy-related protein that has been demonstrated to function as a key autophagy-associated target gene regulated by miR-23b, and suppression of ATG12 significantly increased the radiosensitivity of pancreatic cancer cells, whereas the miR-23b-induced radiosensitivity was eliminated by ATG12 overexpression (125). Studies indicated that ULK served an important role in the autophagy process (126-128). ULK2 and ULK1 are highly homologous, functionally complementary and indispensable for affecting autophagy (129-131). GABARAPL1 is considered as a target of miR-195 and regulates the proliferation, migration, angiogenesis and autophagy of endothelial progenitor cells (132). These three genes (ATG12, GABARAPL1 and ULK2) may serve pivotal roles in the autophagy suppressed by chloroquine diphosphate in pancreatic cancer cells.

To determine relevant miRNAs and elucidate their potential regulatory mechanisms in the autophagy of pancreatic cancer, differentially-expressed miRNAs were analyzed and their target genes were predicted. Subsequently, the results were compared with the target genes detected in the microarray analysis. The co-expressed genes were considered the putative target genes regulated by the abnormally-expressed miRNAs during the autophagy process in pancreatic cancer cells. Bioinformatics analyses were performed on these target genes to identify the signaling pathways and biological processes that were regulated. The results demonstrated that the expression of miR-663a-5p and miR-154-3p was downregulated in the PANC-1 cells treated with chloroquine diphosphate, compared with the control group. miR-663a-5p expression in pancreatic cancer tissues is significantly downregulated and negatively associated with eukaryotic translation elongation factor 1α 2 (eEF1A2) expression due to miR-663 reducing the proliferation and invasion of pancreatic cells by directly targeting eEF1A2 in vitro and in vivo (133). Additionally, miR-663a-5p has been reported to have a low expression in hepatocellular carcinoma (134), non-small cell lung cancer (135) and colorectal cancer (136). Similarly, the present comprehensive analysis indicated that the expression of miR-663a-5p was downregulated in pancreatic ductal adenocarcinoma (PDAC). In the present study, the bioinformatics analyses indicated that the putative autophagy pathways in the pancreatic cancer cells regulated by miR-663a-5p included the aldosterone-regulated sodium reabsorption, Wnt signaling, basal cell carcinoma, fatty acid biosynthesis and PI3K-Akt signaling pathways, of which the Wnt signaling and PI3K-Akt signaling pathways were considered to have a function in the regulation of autophagy (137-143). Currently, whether miR-154-3p is involved in the biological processes of pancreatic cancer cells, including occurrence, development, prognosis and autophagy, remains unknown; however, miR-154-3p was reported to be expressed at a low level in colorectal cancer and was associated with the degree of malignancy of colorectal cancer (144). Similarly, miR-154-3p exhibited low expression in breast cancer and affected the treatment outcome of breast cancer by regulating E2F transcription factor 5 (145). Additionally, the present comprehensive analysis indicated that the expression of miR-154-3p was downregulated in PDAC. Furthermore, the bioinformatics analyses indicated that the putative autophagy pathways in the pancreatic cancer cells regulated by miR-154-3p included the insulin signaling pathway, FoxO signaling, mitogen-activated protein kinase signaling, ErbB signaling, p53 pathway by glucose deprivation and CCKR signaling map pathways. Activation of insulin like growth factor 1 receptor (IGF-1R) signaling antagonizes the decrease in cell viability of human disc cells through the suppression of apoptosis and enhancement of autophagy (146). IGF-1R increases cell viability during hypoxia, which may be dependent on promoting autophagy by suppressing the PI3K/Akt/mTOR signaling pathway (147). The present study demonstrated that a number of genes, including suppressor of cytokine signaling 1 (SOCS1), plasma membrane intrinsic protein 3 (PIP3) and serine/threonine protein phosphatase type 1 α (PP1), which were involved in the insulin signaling pathway, were abnormally expressed when the autophagic level was decreased. Overexpression of SOCS1 suppressed the PI3K/PIP3/Akt signaling pathway with the subsequent PP1 activation, which has been demonstrated to induce autophagy (148,149). Autophagy is inhibited under the condition of oxygen-glucose deprivation (150). Therefore, we hypothesized that miR-154-3p influences the key genes in the insulin signaling or p53 pathways by glucose deprivation, thereby affecting autophagy. Additionally, due to the involvement of CCKR signaling in pancreatic enzyme secretion (151,152), miR-154-3p may serve a key role in the secretion process. Therefore, the data in the present study will provide novel directions for the investigation of pancreatic cancer and its autophagy.

To elucidate which genes inhibited the miRNAs regulating the autophagy of pancreatic cancer cells through competitive binding, potential target associations between circRNA and miRNA, and between lncRNA-miRNA were predicted based on the competitive expression of ceRNA and miRNA. When these results were combined with the differentially-expressed circRNAs and lncRNAs under different autophagic levels detected by the microarrays, numerous ceRNAs that exhibited target associations with miR-663a-5p and miR-154-3p, and negative associations with the expression of the targeted miRNAs under the same changes in the autophagic level were determined. AC024560.2 is an effective lncRNA residing on chromosome 3 that can be used to predict the metastasis of early cervical cancer to lymph nodes (153). In the present study, AC024560.2 was determined to competitively bind to miR-663a-5p and thus regulate the autophagic level of pancreatic cancer cells by inhibiting the expression of this miRNA. For the other ceRNAs, their associations with tumors have not been reported, and further investigations are required.

Since the molecular mechanisms of tumor occurrence and development are complex and the numerous different physiological processes, including differentiation (154) and aging (155), are frequently accompanied by an aberrant autophagic level (156,157), it is not surprising that the interference of autophagy is determined to be one of the therapeutic methods for tumors. The growth of PDAC could be repressed by hydroxychloroquine in vivo and in vitro (158). PDAC was dependent on autophagy and the use of an autophagy inhibitor may become the breakthrough point of treatment (35). However, in the treatment of PDAC with gemcitabine, a high concentration of gemcitabine is required following adding autophagy inhibitors (37). Therefore, autophagy may serve a dual role in the development of PDAC. Investigating the ceRNA network following inhibiting autophagy, of which these genes may serve role in inhibiting autophagy or promoting autophagy, will provide novel ideas for determining new ways to treat pancreatic cancer. On this basis, the biological function and action mechanism of these two miRNAs will be investigated in further studies.

In summary, the miRNAs that were significantly associated with the autophagy of pancreatic cancer included miR-663a-5p and miR-154-3p. These miRNAs were used as the intermediate regulatory sites, upstream to which multiple putative ceRNAs, including AC024560.2, were present and downstream to multiple targeted regulatory genes, including ATG12 and ULK1. The expression of these miRNAs was regulated through ceRNA, miRNA and mRNA interactions, which could serve important regulatory roles in pancreatic cancer autophagy.

In conclusion, autophagy serves an important role in the development and progression of various malignant tumor types, including lung (159), gastric (160) and renal cancer (161). The investigation of differentially-expressed genes in pancreatic cancer cells under different autophagic levels not only clarifies the regulatory mechanism of the autophagic process in pancreatic cancer cells, but also provides novel ideas for the effective diagnosis and treatment of pancreatic cancer. In the present study, microarray technology was employed to confirm that the change in the autophagic level of pancreatic cancer cells was accompanied by various differentially-expressed genes. Following establishing multiple ceRNA, miRNA and mRNA associations, there is more reason to consider that these genes will serve important roles in the diagnosis and treatment of pancreatic cancer, and thus lay the theoretical foundation for subsequent investigations on pancreatic cancer while indicating new research directions. Subsequently, the regulatory functions of these genes at the cellular and tissue levels will be investigated and their specificity and accuracy will be further validated through a series of in vivo and in vitro experiments.