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.

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.
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.
(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.

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

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.
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 6 th 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.
Study design, methods of data analyses and results are shown in Fig 1.

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). miRNAgene 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).

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  Table. These genes were all over-expressed and targeted by the under-expressed micro-RNAs 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

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.

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 immunerelated 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 miR-NAs 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  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. 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.

Innate immune system pathway
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 exocrinelike, 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 overexpression in the pancreatic TCGA dataset. CD53 is targeted by miR-30c, which is underexpressed 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 overexpressed 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 underexpressed 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.
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][62][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.