Serum level of co-expressed hub miRNAs as diagnostic and prognostic biomarkers for pancreatic ductal adenocarcinoma.

Background: Sensitive and specific non-invasive biomarkers are urgently needed in order to improve the survival of patients with pancreatic ductal adenocarcinoma (PDAC), which is the fourth leading cause of cancer-related death. We aim to identify serum hub miRNAs as potential diagnostic and prognostic biomarkers for PDAC.
Methods: A total of 2578 serum miRNA expression data from 88 PDAC patients and 19 healthy subjects were downloaded from the Gene Expression Omnibus database. Weighted gene co-expression network analysis (WGCNA) was constructed and significant modules were extracted from the network by WGCNA R package. Network modules and hub miRNAs closely related to PDAC were identified. The prognostic value of hub miRNAs was assessed by Kaplan-Meier overall survival analysis.
Results: Two modules strongly associated with PDAC were identified by WGCNA, which were labeled as turquoise and brown respectively. Within each module, twenty hub miRNAs were found. At the functional level, turquoise module was mainly associated with tumorigenesis pathways such as P53 and WNT signaling pathway, while the brown module was mostly related to the pathways of cancer such as RNA transport and MAPK signaling pathway. Utilizing overall survival analyses, five "real" miRNAs were able to stratify PDAC patients into low-risk and high-risk groups. Conclusions: The association of specific Hub miRNAs with the development of pancreatic cancer was established by WGCNA analysis. Five miRNAs (mir-16-2-3p, mir-890, mir-3201, mir-602, and mir-877) were identified as potential diagnostic and prognostic biomarkers for PDAC.

Introduction

Pancreatic ductal adenocarcinoma (PDAC) contributes to nearly 80% pancreatic cancer PDAC which is the fourth leading cause of cancer-related death. Compared with the steady improvement in survival for most other cancer types, little advancement have been made in PDAC, which remains a 3% 5-year survival rate 1. By 2030, PDAC will likely surpass colorectal and breast cancers to become the second leading cause of tumor-related death in both men and women 2. Due to the deep anatomical location of the human pancreas, tumor-specific symptoms such as abdominal mass, jaundice, and weight loss typically emerge only after the tumor has reached an advanced stage, and is either unresectable or has already metastasized to the liver or other organs 3. Up to now, CA19-9 is the only clinically available FDA approved serum biomarker for PDAC. However, due to its low sensitivity and specificity, CA19-9 is primarily used for monitoring tumor recurrence after surgery rather than being used as a primary diagnostic and prognostic biomarker 4-7. In order to improve the survival of patients with pancreatic cancer, sensitive and specific non-invasive diagnostic and prognostic biomarkers are urgently needed.

MiRNAs are small (20-24 nucleotides) non-encoding RNAs that regulate gene expression at the transcriptional or post-transcriptional levels 8. The circulating miRNAs bind to targeted messenger RNA (mRNAs) and generally inhibit translation or reduce the stability of those mRNAs 9. MiRNAs could regulate a broad range of biological processes, like cell-cycle, stress response, cell differentiation, apoptosis and invasion. Aberration of expression and functions of miRNAs has been reported to be greatly involved in tumorigenesis and progression 10-12. Recent studies showed that miRNA-143 can promote cell apoptosis and regulate the process of pancreatic cancer by inhibiting the migration and invasion of pancreatic cancer cells. Eight miRNAs (miR-1301, miR-598, miR-1180, miR-155, miR-496, miR-203, miR-193b, miR-135b) were reported as independent predictors for the survival of patients with PDAC 13,14. Therefore, miRNAs are considered to hold the potential as novel diagnostic biomarkers and therapeutic targets of cancers including PDAC.

Weighted gene co-expression network analysis (WGCNA) is an effective modular biological analysis methodology that could help us to identify and screen co-expressed modules and key biomarkers, such as miRNAs and lncRNAs (long non-coding RNAs) 15-18. WGCNA can group miRNAs into a module or network due to their similar expression in association with a given cancer type. In this study, we applied WGCNA methodology to build co-expressed modules and to identify hub miRNAs potentially involved in the pathogenesis and malignant behavior of pancreatic cancer. We attempted to stratify PDAC patients into low-risk and high-risk groups using a specific group of miRNAs identified from this analysis. Data derived from this analysis is supposed to provide statistically convincing evidence that certain hub miRNAs have the potential to become clinically applicable diagnostic and prognostic biomarkers for PDAC.

2 Materials and methods

2.1 MiRNA expression data collection

Microarray dataset of the serum samples of pancreatic cancer patients was downloaded from Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). Dataset GSE85589 was used to construct co-expression networks and identify hub genes in this study. It contained miRNA expression data of serum samples obtained from 19 healthy subjects and 88 PDAC patients. The platform was GPL19117, Affymetrix Multispecies miRNA-4 Array.

The expression data GSE85589 which had been normalized using Robust Multi-Array Average (RMA)19 was downloaded. Then, probe sets were mapped to gene symbols according to the probe annotation files of the GPL10558 platform. Afterwards, the gene expression values were log2 transformed. Reshape2 and ggplot2 were R-based open-source software packages, which conducted a box plot and correlation plot of normalized samples data for pancreatic cancers (PC) and normal subjects (N) (Supplementary Fig ). Then, miRNA probe sets were mapped to gene-symbols according to the probe annotation file of the GPL19117 platform. If multiple probes mapped to a same gene-symbol, the biggest value of all probes that mapped to the same gene symbol was chosen to represent this particular gene.

2.2 Weighted gene co-expression networks and their modules

WGCNA package was a freely accessible R software (R 3.4.3), for the construction of weighted gene co-expression networks 11. Under “WGCNA” environment of R software, by means of goodSamplesGenes function, miRNA data matrix was checked whether they were the good samples and good genes. Then, we applied the pickSoftThreshold function to assess whether the topology of the networks was scale-free, as is required for WGCNA analysis. β was a softthresholding-power parameter that could strengthen strong correlations and penalize weak correlations between genes. Here, the power of β = 1 (scale free R2 =0.89), which was selected to ensure a scale free network (Fig ). Then, a hierarchical clustering tree was constructed with different branches of the tree representing different gene modules 20. Three different ways can be selected to construct the network and identify modules according to different needs. In this study, the one-step function was used for network construction and detection of consensus modules. (Fig )

2.3 Identification of clinical significant modules and hub miRNAs and functional annotation

To identify clinically significant modules, the module eigengene (ME) was calculated for representing and summarizing each module. This measure can be considered to be representative of gene expression profiles and can be used to capture the maximal amount of variation in the module. To quantify module-trait associations, given that we had a summary profile (ME) for each module, we correlated ME with external traits and searched for the most significant associations. This calculation was referred to as the Module-Trait Relationships 21. For inter-modular analysis, we evaluated the Gene Significance (GS) and Module Membership (MM), the latter is also known as eigengene-based connectivity (kME). GS is the absolute value of the correlation between a specific gene and a clinical trait. Using the GS and ME, we can identify key modules that have a high significance for a clinical trait and important module membership 21.

To identify hub miRNAs, the functions of R software module Eigengenes and signedKME were used to identify highly connected intra-modular miRNAs, also called hub miRNAs. In particular, the moduleEigengenes function calculates the eigengene of a module, which is a virtual gene representing the gene expression profile of the entire module. The signedKME function measures the distance of the expression profile of a gene to that of the module eigengene and, thus, quantifies how close a gene is to a module, i.e., it measures the module membership of a gene. Within a module, the genes with the highest module membership scores are considered as hub genes of that module.

2.4 Functional annotation of hub miRNAs

In analogy to previous studies 18, 22, we selected the 20 most connected hub miRNAs for each key module. The miRNet web tool 23 was used to identify the biological pathways, processes, molecular functions and cellular components statistically enriched for the corresponding miRNA target genes. In particular, this tool identifies the enriched KEGG and REACTOME pathways and Gene Ontology (GO) terms based on experimentally confirmed miRNA targets. Since this tool exploits validated miRNA targets, it ensures a higher reliability than tools based on predicted miRNA targets.

2.5 Kaplan-Meier Survival analyses and efficacy evaluation

In order to evaluate the identified hub miRNAs as prognostic PDAC biomarkers, we conducted survival analyses using the SurvMicro web method 24 Survival analyses were performed on an independent PDAC dataset present in The Cancer Genome Atlas (TCGA) (http://cancergenome.nih.gov/) containing miRNA expression and survival data derived from 54 PDAC patients (Supplementary Table 1). The SurvMicro tool divides samples into high-risk and low-risk groups through the median of the prognostic index obtained via a Cox Regression model. After this, it generates the hazard ratios (HR) with relative confidence intervals (CI) and p-values. We aim to find the “real” hub miRNAs associated with the prognosis of PDAC.

3. Results

3.1 Screening the miRNA modules

A total of 2578 miRNAs were collected after data preprocessing. The selected miRNAs in each category were plotted in their corresponding expression heat-maps (Supplementary Fig . WGCNA analysis of the sequencing data from serum miRNA of PDAC patients and healthy subjects revealed a total of five functional modules (i.e., turquoise, brown, blue, yellow and grey module) (Fig ). In line with previous reported approach, the grey module belongs to the special WGCNA category that has no association with the disease being analyzed 12. As shown in the Matrix of Module-Trait Relationships (Fig ), the module represented by a turquoise color showed a positive association (r = 0.59; p = 1e-22) with PDAC, whereas, the brown module had a negative association with PDAC (r = -0.51; p = 3e-08).Using the GS and MM measures, compared with blue and yellow modules, the turquoise and brown modules showed a significantly higher difference between PDAC group and normal group (Turquoise: r = 0.64; p = 4.8e-20; Brown: r = 0.64; p = 1.9e-09). These data justified the selection of the turquoise and brown modules as key miRNA modules for PDAC. (Fig )

3.2 Detection of hub miRNAs and their functional annotations

The signedKME and module Eigengenes functions in “WGCNA” package were applied to determine the highly-connected intramodular miRNAs, i.e hub miRNAs. In particular, the moduleEigengene function calculates the eigengene of one module, which was as a virtual miRNA representing the expression profile of the entire module.

The signedKME function measured the module membership of each miRNA in a special module. The miRNAs with the highest module membership scores are considered as hub miRNAs of that module. Since the turquoise and brown modules were found to be highly associated with the PDAC serum miRNA network, we set out to determine the hub miRNAs that are represented by selecting miRNAs with the highest module membership scores. It is conceivable that these miRNAs may play important roles in the pathogenesis of pancreatic cancer. The top 20 hub miRNAs identified for each module in the PDAC network are listed in Table , and their network of with their experimentally validated target genes are shown in Figure .

In order to functionally annotate the hub miRNAs, we conducted a functional enrichment analysis limited to experimentally confirmed miRNA target genes. By means of the miRNet web tool, the hub miRNAs in each module was analyzed to identify enriched KEGG and REACTOME pathways and Gene Ontology (GO) terms (Table ). According to KEGG and REACTOME and GO function analyses, the turquoise module was mainly enriched in cancer-related terms such as MAPK and WNT signaling, etc, while the brown module was mainly related to the functions including cancer-associated pathways, such as RNA transport.

3.3 Hub miRNAs validation and further distinction

To determine whether the serum candidate miRNAs identified in this study could be used as prognostic biomarkers, we assessed if hub miRNAs in the turquoise and brown modules were connected with the overall survival (OS) of pancreatic cancer patients 18. We performed both multi-miRNA (top 20 hub miRNAs) survival analyses of PDAC patients in each module by using the SurvMicro tool. By applying the top 20 hub miRNAs as input for the SurvMicro tool, we found that the miRNA profiles in all modules and single miRNA successfully stratified PDAC patients into high-risk and low-risk groups as revealed by Kaplan-Meier survival plots (table ). Specifically, the OS time in the high-risk group patients was found at least two-fold shorter compared with those of patients in the low-risk group (HR 6.93 [95% CI 3.03-15.89] p = 4.687e-06 for the turquoise module, HR 3.80 [95% CI 1.94-7.46] p = 9.923e-05 for the brown module) (supplementary Fig ). In addition, we set out to identify “real” miRNAs that were most significantly associated with PDAC survival outcomes by judging whether the P value of each miRNA was less than 0.05. An increased expression of three miRNAs in turquoise module (mir-16-2-3p (HR 2.05 [95% CI 1.08-3.89], p=0.032), mir-890 (HR 2.48 [95% CI 1.30-4.71], p=0.006), and mir-3201 (HR 2.68 [95% CI 1.41-5.09], p=0.003) and two miRNAs in brown module ( mir-602 (HR 2.68 [95% CI 1.41-5.09], p=0.003), and mir-877 (HR 5.58 [95% CI 1.31-23.84], p=0.020) were found to be associated with a significantly better survival (Fig ).

Discussion

Since its introduction in 2005, WGCNA methodology has been proven to be an effective tool for describing the correlation patterns among genes across microarray samples 25. WGCNA is particularly useful in finding modules of highly correlated genes and summarizing many intramodular hub genes that can be extrapolated for the diagnosis and prognosis of cancers. WGCNA involves several sequential steps. Firstly, the absolute values of the connection of paired miRNAs were calculated for the co-expression network. Secondly, an adjacency matrix was built to define the miRNAs strength with each other. A soft threshold parameter was under the weighted network. Topological overlap (TOM) was used to cluster miRNAs into different network modules that combined the adjacency of two direct miRNAs and the connection strengths with other indirect miRNAs. The miRNAs in one module could be characterized with the module eigengene (ME). To find the key modules related to clinical traits, the correlations were calculated between MEs and those clinical traits. The correlation between miRNAs and MEs was applied to distinguish hub miRNAs 21. In this study, we applied this methodology to build co-expressed modules and identify hub miRNAs potentially involved in the tumorigenesis of pancreatic cancer. Among those identified hub miRNAs, we further determined “real” miRNAs that could predict the behavior of PDAC.

Similar work by Siân Jones 26, et al has identified 12 partially overlapping signaling pathways that are genetically altered in the great majority of pancreatic cancer. By means of WGCNA, we found that 20 hub miRNAs were closely associated with, MAPK, P53 and WNT signaling pathways, which are included (WNT and MAPK signaling) or greatly involved (p53 in cell cycle, invasion, DNA damage and cell apoptosis signaling) in the pathways mentioned above. More concretely, target genes of those miRNAs includes TP53BP2 (tumor protein p53 binding protein 2), dkk3 (WNT signaling pathway inhibitor dickkopf-3), and WISP2 (WNT1 inducible signaling pathway protein 2) 27. We also identified another set of 20 hub miRNAs that were closely associated with the functions of PDAC and RNA transport, whose targeted genes include TP53INP2 (tumor protein p53 inducible nuclear protein). Our data suggest miRNAs are key regulators of these core tumorigenic and functional pathways in PDAC.

High serum levels of certain miRNAs might indicate the presence of an underlying malignancy. Mahito Miyamae reported that high level of six circulating oncogenic miRNAs (miR-615-5p, -744, -575, -557, -675, and -550a) could be useful biomarkers for PDAC screening, monitoring recurrence after surgery, and predicting the prognosis of patients with PDAC 28. Schultz et al reported that the area under the curve (AUC) was 0.91 (95% CI, 0.87-0.94) when using whole blood level of 10 miRNAs (miR-26b, miR-34a, miR-122, miR-126, miR-145, miR-150, miR-223, miR-505, miR-636, miR-885-5p) to diagnose PDAC 29. Moreover, the increased miR-16 and miR-196a combined with CA19-9 have been shown to be more effective to distinguish pancreatic cancer from normal individuals 30. Brand et al found that the combination of CA19-9, ICAM-1 and OPG could discriminate PDAC patients from health controls with the 78% sensitivity and 94% specificity 31. In the present study, we identified a list of specific Hub miRNAs that are associated with PDAC tumorigenesis and a list of hub miRNAs that are associated with the development of pancreatic cancer. The association of circulating miRNAs with PDAC tumorigenesis and progression indicates that circulating hub miRNAs could not only be a diagnostic biomarker, but also predict the malignant behavior of the tumor.

To further study the prognostic value of circulating hub miRNAs, we identified 5 “real” miRNAs (mir-16-2-3p, mir-890, mir-3201, mir-602, and has-mir-877) that could stratify PDAC patients into high and low risk groups. miR-15b/16-2 cluster could modulate the expression of cell cycle-regulating genes, CCND2, CyclinD2 and IGF1R, and was deeply involved in proliferation and anti-apoptotic pathways of different kind of malignancies such as leukemia and colorectal cancer 32. Mir-15b was also suggested to promote epithelial-mesenchymal transition (EMT) and metastasis in pancreatic cancers by targeting smurf2 33. MiR-16-2-3p has been reported in several studies as a diagnostic biomarker for hepatocellular carcinoma with a high sensitivity and specificity 34, 35. Mir-890 belongs to a mir-888 genomic cluster located on human X chromosome Xq27.3 that specifically targets RBL1, KLF5, SMAD4, and TIMP2 and functions as pro-oncogenic factors that increased prostate tumor growth in vivo 36, whose function, however, has not been stratified in other cancer types 37. The function of miR-3201 has been still controversial. MiR-3201 expression was significantly reduced in recurrent epithelial ovarian cancer (EOC) 38, indicative of a suppressive role during cancer recurrence. Other study showed high expression of mir-3201 in melanoma patients 39. However, in the current study, serum level of mir-3201, which was suggested to be involved in several of cancer-promoting pathways (table 1), predicts a poor prognosis of PDAC patients. The discrepancies between these studies might arise from the cancer types of different tissue origin, or the sample type, with the former being cancerous tissues and the latter being serum. Mir-602 plays a pro-carcinogenic role in HBV-related hepatocarcinoma by inhibiting RASSF1A 40, and with other three miRNAs, was selected to identify the overall survival of glioblastoma multiforme 41. It can also regulates the expression of sonic hedghog pathway in chondrocytes 42, which was also one of the core signalings during the tumorigenesis of PDAC 26; miR‑877 has been shown to suppress the proliferation and migration of (RCC) cells by modulating the eEF2K/eEF2 signaling pathway 43. The binding site of miR-877-3p was on the promoter site of cancer suppressor gene p16 which alters frequently in tumors. MiR-877-3p targets the p16 promoter and are related to increased risk of a wide range of cancers such as bladder cancer, renal cell carcinoma cells as well as pancreatic ductal adenocarcinoma 44-46. It is conceivable that these 5 miRNAs involved in the tumorigeneis and mailgnant function of PDAC via completely different mechanisms. The conbination of these 5 miRNAs can be a reliable biomarker for predicting the behavior of PDAC.

Table 1

20 Hub miRNAs identified in the PDAC network and Functional annotation of hub miRNAs and enriched KEGG pathways

Module	Top 20 microRNAs	KEGG	Reactame	GO:BP	GO:MF	GO:CC
Turquoise	miR-3201 miR-8084 miR-3128miR-3910 miR-548x-3p miR-1298-3pmiR-32-5p miR-4445-3p miR-5681amiR-890 miR-335-5p miR-603miR-3152-3p miR-16-2-3p miR-3613-3p miR-155-3p miR-548u miR-628-5p miR-130a-3p miR-570-3p	Pancreatic cancer,Pathways in cancer,MAPK signaling pathway,ErbB signaling pathway,Wnt signaling pathway,Prostate cancer,colorectal cancer,Cytokine-cytokine receptor interaction,T cell receptor signaling pathway	Cytokine Signaling in Immune system,Transmembrane transport of small molecules,Axon guidance, Signalling by NGF,Signaling by PDGF,Downstream signal transduction,NGF signalling via TRKA from the plasma membrane,(Signalling by NGF,FGFR1-4,and ERBB2,ERBB4),Developmental Biology,Generic Transcription Pathway	regulation of anatomical structure morphogenesis,regulation of cell proliferation	regulation of transcription DNA-dependent,enzyme binding,zinc ion binding	nucleoplasm
Brown	miR-7b-5p miR-93-5p miR-181b-5pmiR-182-5p miR-106b-5p miR-595miR-602 miR-610 miR-623miR-297 miR-877-5p miR-1225-5pmiR-1469 miR-1470 miR-1914-5pmiR-762 miR-3171 miR-3188miR-4293 miR-4270	Pathways in cancer,Prostate cancer,Insulin signaling pathway,RNA transport,Non-small cell lung cancer,Pancreatic cancer,MAPK signaling pathway,Regulation of actin cytoskeleton	Gene Expression,Cell Cycle,Immune System,Signalling by NGF,Signaling by FGFR1-4,Generic Transcription Pathway,Signaling by the B Cell Receptor (BCR)	cellular macromolecule catabolic process,cell division	nucleotide binding,transition metal ion binding,purine nucleotide binding	cytosol,nuclear lumen

In summary, using WGCNA technology we identified 20 hub miRNAs that are associated with tumorigenesis pathways and 20 hub miRNAs that are associated with functional pathways in tumors. Among the 40 hub miRNAs, 5 “real” miRNAs were able to stratify the PDAC patients into high and low risk groups. Our data indicate certain hub miRNAs can function as diagnostic and prognostic biomarkers for PDAC.

Supplementary figures and tables.

Click here for additional data file.

Table 2

Single hub miRNA survival analysis tested by TCGA PDAC miRNAs data

Module	Hub miRNA	HR	95%CI	P value
Turquoise	miR-16-2-3p	2.06	1.06-3.98	0.032
	miR-3201	2.68	1.41-5.09	0.003
	miR-890	2.48	1.30-4.71	0.006
Brown	miR-877-5p	5.58	1.31-23.84	0.02
	miR-602	2.68	1.41-5.09	0.003