MicroRNA modulated networks of adaptive and innate immune response in pancreatic ductal adenocarcinoma.

Despite progress in treatment strategies, only ~24% of pancreatic ductal adenocarcinoma (PDAC) patients survive >1 year. Our goal was to elucidate deregulated pathways modulated by microRNAs (miRNAs) in PDAC and Vater ampulla (AMP) cancers. Global miRNA expression was identified in 19 PDAC, 6 AMP and 25 paired, histologically normal pancreatic tissues using the GeneChip 4.0 miRNA arrays. Computational approaches were used for miRNA target prediction/identification of miRNA-regulated pathways. Target gene expression was validated in 178 pancreatic cancer and 4 pancreatic normal tissues from The Cancer Genome Atlas (TCGA). 20 miRNAs were significantly deregulated (FC≥2 and p<0.05) (15 down- and 5 up-regulated) in PDAC. miR-216 family (miR-216a-3p, miR-216a-5p, miR-216b-3p and miR-216b-5p) was consistently down-regulated in PDAC. miRNA-modulated pathways are associated with innate and adaptive immune system responses in PDAC. AMP cancers showed 8 down- and 1 up-regulated miRNAs (FDR p<0.05). Most enriched pathways (p<0.01) were RAS and Nerve Growth Factor signaling. PDAC and AMP display different global miRNA expression profiles and miRNA regulated networks/tumorigenesis pathways. The immune response was enriched in PDAC, suggesting the existence of immune checkpoint pathways more relevant to PDAC than AMP.

Introduction

Pancreatic cancer is the fourth most frequent cause of cancer death, worldwide [1,2]. The majority of pancreatic cancers (~96%) comprise exocrine ductal adenocarcinomas (PDAC) [2]. Other pancreatic tumors include periampullary carcinomas; of these, ~12% are adenocarcinomas of Vater ampulla (AMP). AMP patients have a better prognosis (5-year survival of >45%) [3,4] compared to PDAC, mainly due to early disease detection.

The asymptomatic nature of PDAC is often associated with late diagnosis and poor patient prognosis. This disease is characterized by very aggressive and rapid tumor growth and high incidence of distant metastasis [5]. Despite progress in treatment strategies (surgery, chemo and radiation therapies) in the past 5 years, only 24% of patients with PDAC survive more than a year and <5% are expected to survive more than 5 years after diagnosis [6]. More recently developed immunotherapy treatment strategies have failed to succeed for patients with PDAC. Therefore, there is a need to identify clinically useful biomarkers for early disease detection, as well as to further develop combinations of precision therapeutics for patients diagnosed with the late-stage disease [7].

Molecular mechanisms underlying pancreatic oncogenesis involve both genetic and epigenetic changes [8]. The role of microRNAs (miRNAs) has been widely reported as regulatory molecules associated with tumorigenesis [9]. miRNAs are potent gene expression regulators involved in biological processes including embryonic development, differentiation, apoptosis and cell proliferation [10,11]. An increasing number of non-coding RNAs, including miRNAs, have been reported to regulate immune system homeostasis and immune system development and function [12,13]. miRNAs have roles in the regulation of both innate and adaptive immune responses, in which they control early development of immune cell progenitors, maintenance and differentiation and mature immune cell function [12,14]. Considering that miRNAs have been shown as biomarkers for diagnosis, prognosis, and treatment of patients with cancer [15] they may have potential use as combination therapeutic molecules in PDAC, targeting immune system pathways.

Previous studies compared global miRNA expression profiles in PDAC and AMP and showed that miRNA expression profiles are distinct between PDAC and AMP tumors [16-20].

Here we sought to elucidate deregulated mechanisms modulated by miRNAs in PDAC and AMP cancers by experimental analysis followed by application of systems biology approaches. Our study shows that PDAC and AMP display different global miRNA expression profiles, which reflect distinct miRNA-regulated networks and pathways of tumorigenesis. Notably, an enrichment of immune response pathways was found in PDAC, suggesting that immune checkpoint pathways are more relevant to PDAC than AMP. The identification of miRNAs as biomarkers of innate and adaptive immune response in PDAC opens avenues for investigation of novel therapeutic strategies.

Materials and methods

Ethics statement

This study was performed in accordance with the national and international ethical guidelines and recommendations of the Declaration of Helsinki. This study was approved by the Research Ethics Committee of the Faculty of Medicine, UNESP, Botucatu, SP, Brazil (Protocol # 4382/2012-B).

Patient samples

Formalin-fixed, paraffin-embedded (FFPE) tissue samples (paired tumors and histologically normal surgical resection margins) from patients undergoing surgery for resection of PDAC (N = 19) or AMP (N = 6) cancers were retrospectively collected from the Dept. of Pathology, Faculty of Medicine, UNESP, Botucatu, SP, Brazil. Primary pancreatic ductal adenocarcinoma and pancreatic carcinoma of the ampulla of Vater were obtained from patients naive of chemo or radiotherapy before surgery. Exclusion criteria were: patients diagnosed with secondary pancreatic metastasis and with unresectable disease. Patient demographic and histopathological data are shown in Table 1.

Table 1

Demographic and histopathological data of patients with PDAC and AMP.

Variables	PDACN = 19	N (%)	AMPN = 6	N (%)
Age (years)
Median	59		65
Average (range)	59 (38–77)		60 (31 a 73)
Gender
Male	7	36.9	5	83.3
Female	12	63.1	1	16.7
Tumor grade
WD	15	79	4	66.6
MD	2	10.5	1	16.7
PD	2	10.5	1	16.7
Tumor size (T)
T1-T2	7	36.9	2	33.3
T3-T4	12	63.1	4	66.7
Nodal status
Positive (N1)	13	68.5	2	33.3
Negative (N0)	2	10.5	4	66.7
Undetermined (Nx)	4	21	0	0
Distant metastasis
Positive (M1)	0	0	0	0
Negative (M0)	0	0	0	0
Undetermined (Mx)	19	100	6	100
Stage
Ia	1	5.25	2	33.3
Ib	0	0	0	0
IIa	1	5.25	2	33.3
IIb	17	89.5	2	33.3
III	0	0	0	0
IV	0	0	0	0

WD: well differentiated

MD: moderately differentiated

PD: poorly differentiated

Tissue macrodissection and RNA extraction

Formalin-fixed, paraffin embedded (FFPE) tissue sections were obtained from surgically resected tumors (5–10 sections each, 10 μm thick). Tissue sections were prepared for needle macrodissection using the stereomicroscope Leica EZ4 (Leica Microsystems, Wetzlar, Germany), according to a previously reported protocol [21,22]. Tissue macrodissection was performed in order to isolate the total areas of the tumor including stromal cells, separate from adjacent, histologically normal pancreatic epithelia. A representative example of macrodissected areas (tumor and histologically normal pancreatic tissues) is shown in S1 Fig. Total RNA was then isolated from macrodissected samples and purified using the RecoverAll Total Nucleic Acid Isolation Kit (Ambion, Austin, TX, USA), according to the manufacturer’s protocol. This method allows isolation of total RNA, enriching for miRNAs from FFPE specimens. The protease digestion conditions of the RecoverAll kit are designed to release a maximal amount of RNA fragments of all sizes, including miRNAs. RNA quality and quantity were evaluated by spectrophotometry.

Quantification of miRNA expression

GeneChip microRNA 4.0 array (Thermo Fisher Scientific, Waltham, MA, USA) containing probes representing 2,588 miRNAs (miRBase release 21; www.mirbase.org) was assayed according to the manufacturer´s protocol. Briefly, RNA samples were first subjected to poly-A tail incorporation to 3´-end, followed by a second step of ligation reaction and biotinylation of 3´-poly A tail of RNA. Biotin-labeled RNA was hybridized to the GeneChip miRNA 4.0 array cartridge, for each sample, and detected with Avidin-Streptavidin-Phycoerythrin (PE) conjugate, which binds to biotin-labeled RNA with strong affinity. miRNA array cartridges were then placed into hybridization oven trays and trays were loaded into the hybridization oven and incubated at 48°C with 60 rpm rotation for 16 hours. Upon hybridization, each array was filled with an array holding buffer and allowed to reach room temperature before washing and staining. Wash and stain steps were performed using the appropriate fluidics script for cartridge arrays on the fluidics station. Arrays were scanned and data exported for further analysis using the Expression Console software (Affymetrix) for data summarization, normalization, and quality control.

PDAC and AMP data were analyzed independently in order to identify significantly deregulated miRNAs (fold change, FC ≥ 2 and p<0.05) in each tumor subtype compared to their corresponding normal tissues. In order to decrease the small sample size effect of the AMP sample group (N = 6 tumors and 6 paired normal tissues), miRNA analysis for the AMP sample group did not rely in the p-values alone but utilized False Discovery Rate (FDR) p-values [23] to select miRNAs that are potentially relevant in AMP cancers.

Original, raw data are publicly available on Gene Expression Omnibus (GEO), under accession number GSE125179.

miRNA target prediction analysis

Deregulated miRNAs were subjected to target prediction analysis using two different bioinformatic tools: microRNA Data Integration Portal (http://ophid.utoronto.ca/mirDIP/) [24], using as criteria for target selection, results with “very high” (top 1%) and “high” (top 5%) scores for interaction probability. In addition, we included validated interactions in miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/php/index.php) [25].

Validation of miRNA target gene expression in external data sources (TCGA)

Next, predicted miRNA targets were validated using external data retrieved from The Cancer Genome Atlas (TCGA) (https://cancergenome.nih.gov/). Raw RNA sequencing (RNA-Seq) data (read counts) were used to calculate differential gene expression levels in pancreatic carcinoma samples (N = 178) relative to normal pancreatic tissues (N = 4). For AMP subtype, gene expression data were not available in TCGA (date of access: November 6th 2018). We used edgeR package [26] in the default settings for gene expression analysis and generated a list of differentially expressed genes and corresponding fold change (FC) values. Genes were identified as over- or under-expressed with FDR≤0.05.

In addition, expression levels of miRNA target genes were analyzed using unsupervised clustering approaches, according to neoplastic cellularity of PDAC (high purity vs. low purity tumors). Gene expression data were analyzed using Xena Functional Genomics Explorer tool (https://xenabrowser.net/) [27]. Morpheus—Versatile matrix visualization and analysis software (https://software.broadinstitute.org/morpheus/) [28] was used to verify whether PDAC samples would cluster in separate groups according to the expression levels of the 86 genes and epithelial cellularity. Unsupervised hierarchical clustering was applied using Morpheus tool with the following parameters: Euclidean Distance, Average Data, Cluster Genes, and Samples. Significantly deregulated genes in pancreatic carcinoma, which are targeted by the identified miRNAs and with an inverse correlation with miRNA expression (e.g., miRNA under-expression / target gene over-expression) were used in subsequent pathway enrichment and network analyzes, as described below.

Pathway enrichment and network analyses

ToppGene Suite (https://toppgene.cchmc.org/) [29] was used to identify statistically enriched pathways. Protein-protein interaction (PPI) network analysis was performed using Integrated Interactions Database (IID) (http://iid.ophid.utoronto.ca/iid/Search_By_Proteins/) [30] to assembly PPI networks between targets of differentially expressed miRNAs. Using IID, we obtained PPIs among the query genes, selecting for experimentally validated interactions. PPI networks were visualized using NAViGaTOR v.2.3.2 [31], which comprised target genes as the nodes, and PPIs among these genes as the edges.

Study design, methods of data analyses and results are shown in Fig 1.

Fig 1

Flowchart integrating study design and main results.

Results

A subset of miRNAs is down-regulated and modulates validated gene targets within adaptive and innate immune response pathways in PDAC

We identified 20 significantly deregulated miRNAs (FC≥2 and p<0.05) being 15 down-regulated and 5 up-regulated in PDAC compared to matched normal pancreatic tissues from the same patients (Table 2). Interestingly, miR-216 miRNA family (miR-216a-3p. miR-216a-5p, miR-216b-3p and miR-216b-5p) was identified as consistently down-regulated in PDAC. The highest differences were detected for under-expressed miRNAs: miR-217 (FC = -95.45), miR-216b-5p (FC = -62.86) and miR-148a-3p (FC = -33.75). Furthermore, mirDIP analysis with these 20 miRNAs revealed a total of 24,216 interactions filtered for “high” and “very high” interaction probability. Of these, 1,355 genes were validated targets in miRTarBase (S1 Table) and were further mapped on 16 significantly (p<0.05) enriched pathways. The most enriched pathways with a higher number of genes were the innate immune system response (1,312 genes) and adaptive immune response (826 genes) (S2 Fig). Additional pathways involved genes associated with mechanisms of the immune response regulation (S2 Table). miRNA-gene interaction networks of the innate and adaptive immune response are shown in Fig 2A and 2B. miRNA target genes with functions in the adaptive and innate immune system, as identified in our network analysis, were significantly over-expressed in the 178 PDAC samples from the TCGA dataset (Table 3).

Table 2

Deregulated miRNAs in PDAC compared to normal pancreatic tissues.

Down-regulated			Up-regulated
MicroRNA	FC	p-value	microRNA	FC	p-value
hsa-miR-216a-5p	-52	0.000	hsa-miR-210-3p	10.05	0.017
hsa-miR-216a-3p	-5.82	0.000	hsa-miR-708-5p	9.17	0.029
hsa-miR-1827	-2.36	0.000	hsa-miR-106b-3p	4.2	0.031
hsa-miR-216b-3p	-2.14	0.000	hsa-miR-31-5p	12.89	0.04
hsa-miR-5590-3p	-2.09	0.000	hsa-miR-150-5p	7.02	0.047
hsa-miR-217	-95.45	0.001
hsa-miR-216b-5p	-62.86	0.001
hsa-miR-148a-3p	-33.75	0.001
hsa-miR-30c-2-3p	-2.13	0.001
hsa-miR-141-3p	-10.85	0.007
hsa-miR-5195-5p	-2.2	0.008
hsa-miR-130b-3p	-4.82	0.009
hsa-miR-193b-5p	-2.62	0.015
hsa-miR-6736-5p	-2.38	0.023
hsa-miR-3163	-2.17	0.049

FC: fold change. p-value: ANOVA test.

Fig 2

miRNA-gene interaction networks in PDAC for (a) innate and (b) adaptive immune system pathways.

Genes found as significantly over-expressed (FC≥1.5 and p<0.05 and FDR≤0.05) in the pancreatic TCGA dataset are shown in red and labeled by their symbols, and are significantly enriched by under-expressed miRNAs (green). Genes labeled in pink are predicted miRNA targets and interact directly with over-expressed genes.

Table 3

Expression levels of genes involved in innate and/or adaptive immune system pathways, validated in the PDAC—TCGA dataset (N = 178 tumors and 4 normal tissues).

Genes are over-expressed and predicted to be regulated by under-expressed miRNAs.

Innate immune system pathway
Genes	FC	p-value	FDR
CR1	3.54	5.43E-11	5.22E-08
ARHGAP9	2.04	8.09E-06	0.000779105
CD53	1.97	1.08E-05	0.000931297
HVCN1	1.89	3.28E-06	0.000389503
BCL2	1.73	2.87E-05	0.002037914
CYFIP2	1.55	0.000245959	0.011148321
Adaptive immune system pathway
Genes	FC	p-value	FDR
CD22	3.83	1.31E-10	1.13E-07
CD160	3.15	1.35E-23	1.10E-19
CD79A	2.71	0.000257143	0.011590682
PRKCB	2.64	2.34E-06	0.000302879
CD300LG	2.57	0.001662075	0.045301984
SLAMF6	2.14	5.96E-05	0.003631037
PTPRC	2.01	0.000187886	0.009208655
SIGLEC10	1.79	0.000156912	0.008054502
CD226	1.66	0.001735871	0.0470812
HLA-DOA	1.59	0.000543535	0.020331941
Innate and Adaptive immune system pathways
Genes	FC	p-value	FDR
CD36	2.85	1.64E-08	7.18E-06
ITGAL	2.43	4.83E-07	9.37E-05
CD300LB	2.18	4.38E-07	8.69E-05
SELL	2.14	0.000300675	0.013081778
TLR6	1.73	4.17E-05	0.002772864
PIK3CD	1.67	8.46E-06	0.000805096
CYBB	1.67	0.000433473	0.017120097

FC = fold-change. FDR = false discovery rate.

High and low epithelial cellularity PDAC samples tend to cluster in separate groups according to gene expression levels

Our comparative analyses stratifying the TCGA pancreatic tumor dataset into high and low epithelial cellularity tumors included 128 tumors histologically confirmed as pancreatic ductal adenocarcinoma and classified as High or Low cellularity tumors. Of the 128 PDAC tumors, 76 were High and 52 Low Purity/Cellularity samples (data extracted from TCGA). Expression levels of the 86 significantly deregulated genes in High vs. Low Cellularity tumors are shown in S3 Table. These genes were all over-expressed and targeted by the under-expressed microRNAs identified in our study. Unsupervised hierarchical clustering results showed that High and Low cellularity tumors tended to cluster in separate groups, with two enriched, smaller groups containing only High or Low cellularity tumors, and another group containing a mixture of these distinct classifications (Fig 3). In addition, enriched pathways for the 86 genes targeted by miRNAs and deregulated in the PDAC–TCGA dataset are shown in Fig 4.

Fig 3

Heatmap showing unsupervised hierarchical clustering of high and low epithelial cellularity tumors according to gene expression levels (86 over-expressed genes regulated by under-expressed miRNAs).

High and Low cellularity tumors tend to cluster in separate groups.

Fig 4

Enriched pathways for the 86 over-expressed genes targeted by under-expressed miRNAs in high and low epithelial cellularity PDAC samples.

Deregulated miRNAs modulate target gene expression involved in cell growth and differentiation pathways in AMP

We identified 9 miRNAs (8 under- and 1 over-expressed) as significantly differentially expressed (FDR p-value < 0.05) in AMP vs. normal tissues (Table 4). Target prediction analysis using mirDIP for these 9 miRNAs showed 10,787 interactions. Of these, 2,009 were validated interactions for 774 target genes on miRTarBase (S4 Table). These 774 genes were involved in 14 significantly (p≤0.01) enriched pathways (S5 Table). Of these, the most enriched pathways based on the number of involved genes were RAS and Nerve Growth Factor (NGF) signaling. miRNA-mRNA interaction networks for genes in these pathways are shown in Fig 5A and 5B.

Table 4

Deregulated miRNAs in adenocarcinoma of Vater ampulla (AMP).

Under-expressed
miRNA	FC	FDR p-value
hsa-miR-509-5p	-2.27	0.012
hsa-miR-4681	-2.71	0.024
hsa-mir-570	-2.02	0.026
hsa-miR-3925-5p	-2.89	0.030
hsa-miR-20b-3p	-2.19	0.035
hsa-miR-6817-5p	-2.56	0.036
hsa-miR-6766-5p	-2.94	0.036
hsa-miR-3941	-2.28	0.040
Over-expressed
miRNA	FC	FDR p-value
hsa-miR-4429	3.54	0.024

Fig 5

miRNA-gene interaction network in AMP for (a) RAS signaling and (b) Nerve Growth Factor (NGF) signaling pathways. Genes (shown in blue and labeled by their symbols) are significantly enriched by differentially expressed miRNAs (down-regulated: green; up-regulated: red).

Discussion

Recently studies have provided evidence of an immune-related component of PDAC tumors. Here, we provide novel data on the potential role of miRNAs in the regulation of immune-related gene networks including pathways of adaptive and innate immune response involved in PDAC. We validated an inverse correlation in miRNA and target gene expression in these pathways using the pancreatic cancer TCGA dataset; these data suggest a potential role of miRNAs in the regulation of immune-related genes in PDAC tumorigenesis.

Of note, our study used needle macrodissected samples including tumor cells and stromal cells, an important component of the tumor microenvironment. Another study [17] have also identified novel miRNAs associated with PDAC when analyzing tumor and stromal cells, without using any method of enrichment for tumor cells. The analysis of tumor and stromal cells allows the detection of miRNA expression profiles representative of tumor-stroma interactions. This is important, since invasive tumors such as PDAC have a characteristic, extensive desmoplastic stroma, which depends, at large, on regulatory signals to and from tumor cells. The desmoplastic stroma includes inflammation, neovascularization and activation of fibroblasts and cancer-associated fibroblasts, which actively participate in the tumor-stroma crosstalk, driving tumor development and progression [32,33].

Recently, it has been shown that the neoplastic cellularity of PDAC (high vs. low cellularity tumors) impacts the analyses and interpretation of transcriptomic data including coding and non-coding RNAs [34]. Results from our comparative analysis in high vs. low cellularity PDAC samples corroborate the current literature, showing that PDAC tumors are very complex and gene expression levels are associated with neoplastic epithelial cellularity in this carcinoma. Of note, our unsupervised hierarchical clustering results showed higher gene expression levels in low compared to high cellularity tumors. Since low cellularity tumors may contain tumor cells as well as stromal tissue; these data are in agreement with our miRNA data, which are reflective of miRNA alterations in tumor cells and stromal cells in PDAC.

Overall, our study identified a smaller number of deregulated miRNAs in PDAC, compared to one of the largest studies in PDAC [17]. Scientific and technical differences may explain these findings, such as the different number of patients, distinct platforms and methods of data analyses, and the nature of control samples used in both studies. Nevertheless, and more importantly, when we compared our data with the data by Schultz et al. [17], we observed that 6 miRNAs (miR-130b-3p, miR-141-3p, miR-148a, miR-216a, miR-216b, and miR-217) were commonly under-expressed in both studies. Considering that patients in both studies were from different geographical locations, we highlight this commonly altered subset of miRNAs as relevant for PDAC.

Recently, integrative transcriptional profiling analysis of primary PDAC, pancreatic human and mouse cell lines defined three PDAC subtypes: classical, quasimesenchymal and exocrine-like, which show differences in outcome and therapeutic responses [35]. In another study [36], basal-like and classical pancreatic tumor subtypes were identified. Additionally, an immunogenic subtype of pancreatic cancer has been described [37], in which tumors harbor upregulated immune networks and distinct molecular pathways that may explain the evolution of this different disease subtype. An immune stromal component was validated in PDAC, containing a high expression of CD37, CD53, CD4, and CSF1R [38]. We validated CD53 over-expression in the pancreatic TCGA dataset. CD53 is targeted by miR-30c, which is under-expressed in our PDAC dataset.

PDAC has a characteristic complexity of inflammatory and regulatory immune cells in tumor cell infiltrates in crosstalk with the tumor microenvironment, which is composed of a variety of immune cells [39,40]. In the tumor microenvironment, cytokine production can promote cancer cell growth [41]. We found a number of cytokines regulated by miRNAs in innate immunity, including IL-6, which is a direct target of miRNAs miR-216-5p and miR-217, CXCL8/IL-8, targeted by miR-5590-3p, IL-11 regulated by miR-1827. Under-expression of the miR-216 cluster (miR-216a-3p, miR-216a-5p, miR-216b-3p and miR-216b-5p) was consistently observed in our PDAC samples compared to normal pancreatic tissues. Significant reduction in expression of the miR-216 cluster has been reported in pancreatic cancer and induced expression of these miRNAs decreased tumor cell aggressiveness, indicating that these miRNAs have a tumor suppressor effect [42].

In pancreatic cancer, high expression levels of IL-6, a cytokine produced by a variety of cells including immune and pancreatic tumor cells [43], and IL-1β were correlated with poor overall and progression-free survival in pancreatic cancer patients [44]. IL-6 gene is strongly over-expressed in pancreatic cancer compared to normal tissues [45,46] and IL-6 protein has strong cytoplasm expression in pancreatic tumor cells [47]. IL-6 is also highly expressed in serum from PDAC patients compared to non-diseased individuals and associated with disease aggressiveness and inflammatory response [48] and pancreatic cancer-related cachexia [47]. In our network analysis, we show miRNA-gene interactions highlighting known mechanisms such as the signal transduction activated by IL-6 that in turn activate JAK2 and STAT3. Importantly, a recent study demonstrated that silencing of IL-6 leads to increased sensitivity to gemcitabine treatment [49].

IL-1β, IL-4, and IL-8 interact with other miRNA-regulated targets within the innate immune response network. These cytokines in addition to VEGF, TGF, and IL-10 were consistently over-expressed in different sources of samples from PDAC patients [50]. CXCL-12 production is regulated by a number of miRNAs found under-expressed in our PDAC data: miR-30c-2-3p, miR-148a-3p, miR-216b-5p and miR-1827. CXCR4 was predicted to have a strong interaction with JAK3, which is regulated by miR-30c-2-3p and miR-217, which are under-expressed in PDAC.

The adaptive immune response pathway was correlated with a number of down-regulated miRNAs in our PDAC dataset: 8/20 (40%) miRNAs identified in our study directly target genes in this pathway. Specifically, members of miR-216 family (miR-216a-5p, miR-216a-3p and miR-216b-5p) and miR-217 target CD antigens including LILRB1 and LILRB4 [51], CD36, CD160 [52], CD226 [53], CD300LB and SIGLEC10 [54] as well as toll-like receptor TLR6. In the context of the tumor microenvironment, gene expression alterations may be dependent on post-transcriptional regulatory mechanisms such as miRNAs, which have been shown to regulate the recruitment and activation of immune cells to the tumor [12].

Considering that the PDAC tissues contain a variety of cells including tumor cells, inflammatory cells, fibroblasts, we searched for the specific cellular localization of miRNAs. For this, we used the Functional Annotation of Mammalian Genome (FANTOM5) consortium data. This consortium provides small RNA sequencing data on 396 human samples across 118 cell types obtained from different organs [55]. We searched this database for all miRNAs found as significantly under-expressed in tumor vs. normal (Table 2). The results showed enriched Gene Ontologies for a subset of miRNAs: miR-216a-3p, miR-216b-5p, miR-217 and miR-130b-3p are mainly expressed in endothelial cells of vascular tree; miR-141-3p and miR-3163 are mainly expressed by epithelial cells, and miR-148a is mainly expressed by mast cells, myeloid leukocytes, and secretory cells. Of note, mast cells are present in the tumor stroma, contributing to the inflammatory microenvironment that impacts tumor behavior [56,57]. In the context of our findings, mir-148a under-expression in tumor and stroma potentially leads to increased expression of its target genes such as BCL2, GRIN2B, and SELL, which play roles including adaptive immune system regulation. Recent studies showed that tumor cells are capable of evading immune response by increasing the expression of genes that encode inhibitory molecules or adhesion molecules. Interestingly, SELL over-expression has been associated with differential immune cell infiltration in lung cancer [58]. In pancreatic ductal adenocarcinoma, ITGAL, which encodes another adhesion molecule and is also shown in our miRNA-mRNA networks, was over-expressed in tumor cells. Interestingly, ITGAL has been associated with enriched pathways of granulocyte adhesion and diapedesis and leukocyte extravasation signaling [59]. In conclusion, the identified miRNAs may play important roles in tumor cells modulating important oncogenic pathways being reflective of tumor-microenvironment interactions, particularly by modulating immune tumor response.

In AMP cancers, we highlight RAS, NGF and EGFR and axon guidance signaling pathways, since identified miRNAs modulate the expression of target genes with overlapping roles in all three pathways. Altered expression of genes in these pathways is associated with deregulated expression of 9 miRNAs: over-expression of miR-4429 and under-expression of miRNAs miR-20b-3p,miR-509-5p, miR-570-3p, miR-3925-5p, miR-3941, miR-4681, miR-6766-5p, and miR-6817-5p.

Interestingly, miR-323a-3p was also reported as over-expressed in pancreatic ampullary tumors and found to participate in interaction networks, modulating gene expression including KRAS, with roles in axon guidance pathway [60]. Results from this and our study suggest alternative miRNA regulatory networks modulating similar pathways in AMP cancers. Additionally, the NGF pathway has been previously reported to be induced by inflammatory cytokines in experimental models of liver disease [61-63]. NGF has been shown to alter the expression of several genes implicated in ovarian oncogenesis and a potential mechanism of NGF-regulated protein synthesis is by post-transcriptional miRNA regulation [64]. Notably, axon guidance, NGF and EGFR signaling mechanisms share a number of genes that have been already demonstrated to play important roles in tumor development and progression, such as the regulatory subunits of the phosphoinositide 3-kinases, PI3Ks, a family of lipid kinases that coordinate several cellular functions including proliferation and survival and are deregulated in different cancer types, including pancreatic cancer [65,66]. Genes within the EGFR signaling pathway have been shown to be frequently altered at both genomic and protein levels, contributing to tumor growth and progression of ampullary pancreatic carcinoma [67]. Axon guidance has been reported in the regulation of numerous biological processes including cancer [54,68]. Genes regulating axon guidance were frequently mutated and showed copy number alterations in PDAC at early stages [69]. The axon guidance factor netrin-1 has been shown to mediate effects of MUC4 gene in neural invasion in PDAC [70]. Considering that our subset of patients with AMP had early-stage tumors, these data suggest that axon guidance pathways may be involved in early tumor development in different histological subtypes of pancreatic carcinoma.

In conclusion, our study identified miRNA changes in tumor and stromal cells of PDAC and AMP; such changes likely contribute to modulate the translational efficiency of corresponding target genes; however, we highlight that there are other mechanisms that interplay for translational regulation influencing gene expression. Indeed, translational control of gene expression occurs through a highly complex and dynamic interplay between mRNA sequences (cis-regulators) and the translational machinery, or canonical translational factors [71]. In addition to these factors, miRNAs are among the trans-regulator molecules that interact with the translational machinery core, to modulate translation of mRNA [Reviewed in [72]].

In addition, while functional studies may lead to a deeper understanding on the role of identified miRNAs on specific pathways, the coordinated functions of multiple miRNAs and their targets into regulatory networks have yet to be elucidated in the different subtypes of pancreatic cancer.

Our results suggest that miRNA-modulated pathways of the adaptive and innate immune response may represent an opportunity for improved therapeutics to benefit patients with pancreatic ductal adenocarcinoma.

Conclusions

Our data provide a comprehensive picture of distinct pathways modulated by miRNAs in different subtypes of pancreatic tumors. Adaptive and innate immune response pathways modulated by miRNAs may be useful as potential therapeutic strategies specifically benefiting patients with PDAC.

Representative example of macrodissected areas: Hematoxylin and Eosin (H&E) stained section showing in (A) adjacent normal tissue; normal pancreatic ducts are shown by the arrows and (B) pancreatic ductal adenocarcinoma cells and stromal component. Magnification: 200 μm.

(TIF)

Click here for additional data file.

Pathways with a higher number of miRNA-target genes in PDAC were the innate and adaptive immune responses.

Graphical representation of miRNA target genes (Y-axis) according to significantly enriched pathways (p<0.05) (X-axis).

(TIF)

Click here for additional data file.

Pathways with a higher number of miRNA-target genes in AMP tumors were RAS and Nerve Growth Factor (NGF) signaling.

Graphical representation of miRNA target genes (Y-axis) according to significantly enriched pathways (p<0.01) (X-axis).

(TIF)

Click here for additional data file.

Predicted miRNA target genes in PDAC, validated in miRTarBase.

(XLSX)

Click here for additional data file.

miRNA modulated pathways in PDAC.

(XLSX)

Click here for additional data file.

Expression levels of 86 deregulated genes in PDAC stratifying by epithelial cellularity.

(XLSX)

Click here for additional data file.

Predicted miRNA target genes in AMP, validated in miRTarBase.

(XLSX)

Click here for additional data file.

miRNA modulated pathways in AMP.

(XLSX)

Click here for additional data file.