Next-generation sequencing reveals novel differentially regulated mRNAs, lncRNAs, miRNAs, sdRNAs and a piRNA in pancreatic cancer

Background Previous studies identified microRNAs (miRNAs) and messenger RNAs with significantly different expression between normal pancreas and pancreatic cancer (PDAC) tissues. Due to technological limitations of microarrays and real-time PCR systems these studies focused on a fixed set of targets. Expression of other RNA classes such as long intergenic non-coding RNAs or sno-derived RNAs has rarely been examined in pancreatic cancer. Here, we analysed the coding and non-coding transcriptome of six PDAC and five control tissues using next-generation sequencing. Results Besides the confirmation of several deregulated mRNAs and miRNAs, miRNAs without previous implication in PDAC were detected: miR-802, miR-2114 or miR-561. SnoRNA-derived RNAs (e.g. sno-HBII-296B) and piR-017061, a piwi-interacting RNA, were found to be differentially expressed between PDAC and control tissues. In silico target analysis of miR-802 revealed potential binding sites in the 3′ UTR of TCF4, encoding a transcription factor that controls Wnt signalling genes. Overexpression of miR-802 in MiaPaCa pancreatic cancer cells reduced TCF4 protein levels. Using Massive Analysis of cDNA Ends (MACE) we identified differential expression of 43 lincRNAs, long intergenic non-coding RNAs, e.g. LINC00261 and LINC00152 as well as several natural antisense transcripts like HNF1A-AS1 and AFAP1-AS1. Differential expression was confirmed by qPCR on the mRNA/miRNA/lincRNA level and by immunohistochemistry on the protein level. Conclusions Here, we report a novel lncRNA, sncRNA and mRNA signature of PDAC. In silico prediction of ncRNA targets allowed for assigning potential functions to differentially regulated RNAs. Electronic supplementary material The online version of this article (doi:10.1186/s12943-015-0358-5) contains supplementary material, which is available to authorized users.


Background
A five year survival of around 4% makes pancreatic ductal adenocarcinoma (PDAC) the fourth leading cause of cancer-related deaths worldwide [1]. To better understand the aggressive growth and the poor response of PDAC to chemotherapeutic agents, studies are required that focus on molecular mechanisms underlying pancreatic cancer development and progression [2].
Somatic mutations, alterations in coding-and noncoding RNA expression as well as in the methylome of pancreatic cancer cells have been intensively studied [3]. The application of NGS technologies such as RNAseq, whole exome sequencing or bisulfite sequencing to PDAC samples provided an unbiased view on genetic and epigenetic alterations [4][5][6]. However, studies investigating the small ncRNAome of PDAC tissues and healthy pancreas utilizing NGS methods are rare.
Studies utilizing microarrays revealed that deregulated microRNAs (miRNAs) have an impact on coding-gene expression in PDAC [7][8][9][10][11][12]. MiRNAs are a class of small non-coding RNAs (sncRNAs) that can repress gene expression [13]. Due to the known technical limitations of microarrays or real-time PCR most studies were focused on a fixed set of miRNA targets and many other sncRNA types have not been implicated in PDAC, but their deregulation and contribution to cancer progression has been described for other cancer types [14,15]. These sncRNAs include small nucleolar-derived RNAs (sno-derived RNAs, sdRNAs), functioning like miRNAs, or regulating splicing and translation [16], as well as piwi-interacting RNAs (piRNAs), that are associated with chromatin organization, messenger RNA stability and genome structure [17].
Similarly, only few studies have reported altered expression of long non-coding RNAs (lncRNAs) in PDAC [18]. LncRNAs represent a diverse class of modestly conserved, polyadenylated, non-protein-coding RNAs with essential roles in tumorigenesis [19]. LncRNAs comprise, among others, long intergenic non-coding RNAs (lincRNAs) and natural antisense transcripts (NATs). NATs, which are transcribed from the opposite ("anti-sense") strand of a protein-coding gene can either stabilize or destabilize the expression of their sense partner [20], lincRNAs act as competitive endogenous RNAs (ceRNAs) and sequester miRNAs ("miRNA sponge") [21], target chromatin modification complexes or RNA-binding proteins to alter gene expressing programs [22].
Here, we used high-throughput NGS-based technologies, namely Massive Analysis of cDNA Ends (MACE) and small RNA-sequencing (sRNA-seq), to characterize the complete coding-and non-coding transcriptomes of tissues from six PDAC patients and five normal controls. We merged the expression pattern obtained by the two techniques to gain insights into the altered miRNA regulation of coding gene expression in PDAC as compared to normal pancreatic tissues. Additionally, we provide evidence for differential expression of a piRNA and several sdRNAs, lincRNAs and NATs in PDAC.
Overall, 396,542,460 sequences were generated. The amount of high-quality reads used for the analysis ranged from 4,150,706 reads in the MACE library for cancer patient P3 to 11,007,400 in the sRNA-seq normal pancreas library N2 (see Table 1). Robust expression was detected for 13,614 coding genes and 432 lncRNAs (Additional file 1: Table S1) whereas 1,961 mRNAs and 43 lncRNAs of these were significantly differentially expressed between cancer and control tissues (Additional file 2: Table S2). Unsupervised hierarchical clustering, principal component analysis (PCA) and Pearson's moment product correlation coefficients (PCCs) revealed a clear separation between the groups based on the sequencing results ( Figure 1A-C). Notably, 70% of all reads across control libraries mapped to 25 genes (PNLIPRP1, CELA3B, CPA2, CTRL, GP2, CPB1, CTRC, RBPJL, KLK1, PLA2G1B, CELA2A, CEL, GNMT, CELA3A, CPA1, PRSS3P2, PRSS3, CLPS, PNLIP, SLC39A5, SPINK1, CTRB1, TMED11P, PRSS1, GATM), encoding pancreatic acinar cell secretoryand related proteins ( Table 2). In cancer libraries, only 3% of reads were annotated to these genes. Their homogenous, robust expression across all controls and strong downregulation in all PDACs underlines the existence of normal pancreas function in control libraries, which is lost in PDAC tissues. The significantly up-/downregulated lncRNAs, as determined by MACE, are listed in Table 3.
From 921 sncRNAs with robust expression (Additional file 3: Table S3), 123 were significantly differentially expressed between the two conditions (Additional file 4: Table S4). Similar to MACE, sRNA-seq allowed a clear separation of cancer and control tissues by PCA, PCCs, and unsupervised clustering (Additional file 5: Figure S1). The twenty most significantly up-/downregulated sncRNAs are listed in Table 4. Among all differentially regulated sncRNAs, we found 104 mature miRNAs, 18 sdRNAs and one piRNA. An overview of differential RNA expression across the genome is given in Figure 2.
Enrichment analysis of protein-coding genes for identification of biological and functional differences Gene ontology (GO)-enrichment analysis of downregulated genes in PDAC revealed a total of 39 significantly enriched GO terms (Additional file 6: Table S5) across all three GO categories. Most of these terms were linked to normal pancreas function, as e.g. "Digestion" (FDR = 1E-5), "Steroid metabolic process" (0.002) or "Triacylglycerol lipase activity" (0.014). In addition, other GO terms corresponding to translation of mRNAs into proteins, as e.g. "Translational elongation" (5.6E-15), "Cytosolic ribosome" (3.2E-8) or "Ribosomal subunit" (2E-5) were enriched, indicating a loss of normal pancreas function.
Upregulated genes in PDAC were enriched in 208 GO terms (Additional file 7: Table S6), including "Immune response" (1.2E-17), "Cell proliferation" (4.45E-5) and "Cell migration" (9E-5). The second set represents an enhanced proliferative potential, the third a high metastatic potential of the cancer cells. The most significantly enriched GO term was "Extracellular matrix" (2.2E-30), corresponding to genes involved in fibrogenesis, such as collagens and fibronectin, as well as TGFB and genes related to this pathway. Several categories were related to a sustained angiogenesis, like "Vasculature development" (2E-8) or "Blood vessel development" (8E-8). Taken together, GO-analysis confirms a loss of normal pancreatic function in the tumor tissues, in combination with increased proliferative potential, extracellular matrix (ECM) activation, blood vessel formation, metastatic spread and the potential to escape the immune system.
A total of 78 sncRNAs (74 miRNAs, 4 sdRNAs) showed significant upregulation in PDAC. Several of these were previously implicated in pancreatic cancer development (e.g. miR-21, 143, 222, 155, 10a, 874) others have not been For each sample of control (N) and PDAC (P) tissues the number of raw sequencing reads, PCR duplicate-free reads, reads mapped to regions of interest (exons -MACE or small ncRNAs -sRNA-seq) and the number of deregulated RNAs is provided for MACE and sRNA-seq libraries, respectively.
shown to be upregulated in PDAC before (e.g. miR-31, 511, 2355). The expression of all differential miRNAs is visualized in Additional file 8: Figure S2.

Overexpression of miR-802 decreased TCF4 protein expression
Considering the absent expression of miR-802 in PDAC tissues and the in silico predicted binding sites of TCF4, we re-expressed miR-802 in MiaPaCa PDAC cells and assessed TCF4 expression. First, we induced miR-802 reexpression in MiaPaCa cells transfected with PCMV-MIR-802 ( Figure 4a). Highly elevated levels of miR-802 were observed 24 h after transfection. To assess the effect of miR-802 on TCF4 protein levels, we harvested transfected MiaPaCa cells and analysed the proteins by western blot analysis ( Figure 4b). Here, TCF4 protein levels decreased to 67% as compared with samples transfected with the negative control ( Figure 4c).

Validation of ZEB1 expression at the protein level
Since we hypothesize that miR-802 regulates ZEB1 expression via TCF4, we analysed expression levels of ZEB1 by immunohistochemistry in samples of human PDAC (n = 10) and normal pancreatic tissues (n = 10). In normal pancreatic tissue, ZEB1 was sparsely seen in periacinar cells (e.g. stellate cells). As depicted in Figure 5, we detected ZEB1 in PDAC samples in stromal cells within desmoplastic areas, but epithelial tumor cells did not express ZEB1. In accordance with previous observations [29] we detected ZEB1 in all analyzed pancreatic cancers only in the stromal compartment but not in epithelial cells. This observation emphasizes ZEB1 as a mesenchymal differentiation marker. We speculate that the newly identified miRNA802 might be involved in the regulation of ZEB1 and consequently might promote the mesenchymal character of pancreatic cancer.

MiRNA-mRNA interaction network analysis of differentially expressed genes
To test whether miRNAs are involved in tumor-specific functional categories detected by GO-enrichment analysis, we exemplarily created a miRNA-mRNA interaction network for the term "Cell motion" (Figure 6). The network contains interactions between gene products of upregulated transcripts in PDAC, as well as downregulated miRNAs whose loss might cause the upregulation of their target genes, predicted by at least three independent miRNA-mRNA interaction tools.

Confirmation of differentially expressed genes, lincRNAs and miRNAs by qPCR
To validate sequencing results obtained by sRNA-seq and MACE, the expression of seven candidate genes, upregulated in PDAC (TCF4, ZEB1, CD1A, FOXL1, GPR87, KLK7, CTHRC1), one down-and three upregulated miR-NAs (miR-802, 135b-5p, 145-5p, 103a-3p) as well as two The normalized logarithmic expression of the ten most highly expressed genes in normal pancreas tissues for all five control-(N) and six PDAC (P) sequencing libraries. All genes encode pancreas-specific proteins and are homogenously, highly expressed across all control libraries and more than~100-fold decreased in all PDAC tissues.  Figure 7.
Except for miR-103a-3p, all tested RNAs were below a significance threshold of p < 0.05 (Wilcoxon's rank sum test) when comparing expression of control and tumor samples. This is consistent with sequencing results, where similarly only miR-103a-3p did not reach the level of significance.
CD1A, FOXL1, GPR87 and KLK7, mRNA expression was undetectable by qPCR in normal pancreatic tissue. This is consistent with sequencing results, where no or at most two reads (normalized) were annotated to these transcripts. Similar to sequencing, all four genes were robustly expressed in PDAC tissues, with C t values ranging from 24.8 to 31.1 and a median C t of 27.3.
Consistent with MACE results, CTHRC1, TCF4 and ZEB1 were significantly higher expressed in PDAC compared to normal pancreas tissues, with median normalized C t (ΔC t ) values between −1.8 and 1.8 in cancer-and 3.1-5.1 in normal tissues. The expression differences (log2 ratio normal/tumor tissue) measured by qPCR (negative ΔΔC t ) highly agree with the log2 fold-changes of the sequencing results (PCC = 0.86).
Similarly, downregulation of miR-802 was validated by qPCR, with a median negative normal/tumor ΔΔC t of 8.2 and p = 0 (Wilcoxon's rank sum test), which is consistent with the log2fc and p-value estimated by sRNAseq (10.9, 0). Similar agreements between sequencing and qPCR results were obtained for the other three miRNAs.

External validation of MACE data by microarray
Badea and colleagues [30] investigated 36 whole tumorand adjacent normal pancreatic tissue samples by codinggene microarray analysis. We compared the 53 most significantly upregulated genes in tumors from their publication with the MACE sequencing results based on logarithmic fold-change (Additional file 11: Figure S4) and statistical significance (p-value). Only five of the genes (DCBLD1, PGM2L1, PDGFC, COX7A1, LY6E), were not significantly upregulated (p < 0.05) in the MACE data, whereas 48 genes showed a significant upregulation as detected by both methods. The PCC between log2fc of MACE and microarray data (0.61) indicates a strong correlation between the results and underlines the reliability of the approach.

Discussion
Our study investigated the coding-and non-coding transcriptomes of six PDAC patients and five healthy pancreatic control tissues. We detected 1,961 mRNAs, 43 lncRNAs and 123 sncRNAs as differentially expressed between the groups. Among these are several coding and non-coding RNAs without previous implication in pancreatic cancer development, most prominently miR-802 which is strongly downregulated in PDAC. Bioinformatic and functional analysis revealed post-transcriptional regulation of TCF4 protein levels by miR-802. Differential regulation of four miRNAS (miR-802, miR-135b-5p, miR-145-5p, 103a-3p), seven genes (CD1A, FOXL1, GPR87, KLK7, CTHRC1, TCF4, ZEB1) and two lincRNAs was confirmed by qPCR.
MiR-802 is the third most significantly repressed miRNA in PDAC, besides the tumor suppressor miRNAs miR-216 and miR-217 that -among others -target KRAS, PTEN, and SMAD7 [29]. MiR-802 is mainly expressed in pancreatic acinar cells [31], which may be the cells of origin for pancreatic preneoplastic lesions and pancreatic cancer [32]. A significant downregulation of miR-802 is observed in mice with ethanol-induced chronic pancreatitis, which predisposes to pancreatic cancer [33]. In contrast, sRNA-seq of pancreatic cyst fluids from low-grade benign and high-grade invasive lesions revealed thirteen enriched miRNAs, among these miR-216, miR-217, and miR-802, in the cyst fluids derived from invasive carcinomas [29]. The reason for the inverse correlation between the expression levels of these tumor suppressor miRNAs in body fluids and tumors currently remains unexplained.
Since no previous studies have reported downregulation of miR-802 in pancreatic cancer, validated targets are rare. Nevertheless, miR-802 targets were identified in two other cancer-types: osteosarcoma and lung cancer [34,35]. In contrast to PDAC, miR-802 is upregulated in both cancers. MiR-802 elevation promotes proliferation of lung carcinoma cell lines by targeting the tumor suppressor gene MEN1. Similarly, cell proliferation was promoted by   miR-802 in osteosarcoma, where the gene encoding p27, a negative cell-cycle regulator, is a direct target. In hepatocellular carcinoma miR-802 is more than 100-fold downregulated, but no targets have yet been identified [36]. Bioinformatic in silico prediction points to Wnt signalling related transcription factor TCF4 mRNA (~12-fold upregulated in PDAC, FDR = 5E-7) as a direct target of miR-802. To validate the generated in silico predictions, we re-expressed miR-802 in the PDAC cell line MiaPaCa and analysed TCF4 protein expression. After re-expression of miR-802, we observed a 30% reduction of TCF4, indicating a direct impact of miR-802 on TCF4 regulation.
Furthermore, the regulation of miR-181a/b expression has been associated with TCF4 expression in hepatocellular carcinoma [38]. Consistent with these observations, our study validates the upregulation (4.6-79-fold, FDR 0.002-4.9E-7) of these six miRNAs with important functions in PDAC development [39,40] that have TCF4 binding sites in their promoter.
To test whether specific miRNAs contribute to the metastatic potential of pancreatic cancers, we used omiRas to predict interactions between downregulated miRNAs and upregulated genes from the GO term "Cell motion". The miRNA with the highest number of targets was miR-130b. Previous studies confirmed downregulation of miR-130b in pancreatic cancer and functional tests revealed that overexpression of miR-130b remarkably inhibited the invasiveness of pancreatic cancer cells [43]. MiR-130b loss (8-fold downregulated in PDAC, FDR = 8.2E-7) might lead to the upregulation of metastasis associated key oncogenes MET, NRP2, ITGA4, THBS1 and SPOCK1 [44,45]. Recently published data describe the impact of miR-130b in metastasis formation [46] and therefore validates the approach of in silico prediction of miRNAs.
In contrast to miRNAs, lncRNAs have just recently moved into the focus of cancer research [47,48]. Nevertheless, studies examining the role of lncRNAs in specific oncogenic processes are limited to date.
Liu and colleagues reported an increased expression of the lncRNA MALAT1 in PDAC compared to adjacent normal pancreatic tissues [49]. Furthermore, they described a significant correlation between MALAT1 expression levels and tumor size, tumor stage, invasion, and disease-free survival [49]. In our analysis, however, MALAT1 upregulation in PDAC was not found. This is further supported by the pancreatic expression database (PED), where MALAT1 was also not reported as differentially expressed [50]. B) Co-expression analysis of significantly upregulated transcription factors that harbour predicted miRNA binding sites in their 3′ UTRs for one of the ten most significantly upregulated miRNAs (miRNAs that have no seed sequence for a TF UTR not shown). Significant (p < 0.01) correlations are indicated by a dot, positive correlations are marked in blue, negative correlations in red. The more significant the correlation, the larger the dot size. Sequence complementarity between an UTR and a miRNA is indicated by an "S". C) Expression of TCF4, ZEB1 and miR-21 across all control (C) and PDAC (P) samples (significant positive correlation) as well as the expression of miR-802 (significantly inversely correlated). The normalized expression for each gene/miRNA is given in log2-scale for each sample.
Similar to MALAT1, an increased expression of HOX antisense intergenic RNA HOTAIR has been found in pancreatic tumors [51]. According to the PED, no significant upregulation of HOTAIR is found in PDAC tissues, which is consistent with our results. Overexpression of PVT1 in the pancreatic cancer cell line ASPC-1 resulted in decreased gemcitabine sensitivity [52]. In this regard, we found an approximately 6-fold upregulation of lncRNA PVT1 in pancreatic cancer.
Besides PVT1 we provide evidence for a deregulation of other 42 lncRNAs in PDAC. Among these LINC00152 is, like in PDAC, overexpressed in gastric cancer (GCC), and its high expression correlates with increased invasion [53]. Xia and colleagues speculate about LINC00152 functioning as a competing endogenous RNA (ceRNA) that sequesters miR-18a-5p, 195-5p, 139-5p and miR-31-5p in GCC and thereby influences i.a. THBS1 expression [21]. We report LINC00152 overexpression in pancreatic cancer (5-fold, FDR = 7.1E-7), whereas all miRNAs with binding sites in the transcript are not significantly deregulated compared to control tissues. This indicates that increased LINC00152 expression might decrease the availability of particular miRNAs by competing for their binding and consequently lead to an upregulation of the miRNA target genes, such as THBS1 (5-fold upregulated in PDAC, FDR = 7.7E-6).
Overexpression of snoRNAs is a common feature in breast and prostate cancer [15], and contributes to tumorigenicity in vitro and in vivo. SnoRNA U50 is downregulated in prostate cancer and potentially functions as a tumor suppressor in other cancer types [54]. Several reports describe that snoRNAs are further processed to generate smaller fragments (sdRNAs) with miRNA-like functionality [55]. Currently, there is no evidence for snoRNAs/sdRNAs involved in pancreatic cancer development. To our best knowledge, this is the first report of differential regulation of snoRNAs/sdRNAs in PDAC. The most significantly regulated sdRNA (34 bps long) is from sno-HBII-296B (SNORD91B), which is approximately 5fold downregulated in PDAC (FDR = 5.2E-5). However, its functional role and that of other differentially expressed snoRNAs/sdRNAs remains unclear; similarly, the functions of piRNAs has not been fully elucidated. PiRNAs are usually involved in germline development, silencing of  transposons, and maintenance of DNA integrity [56]. Upregulated expression of piR-651 has been described in several cancer cell lines, where it promotes cell growth and might serve as a marker for cancer diagnosis [56]. Here we report the downregulation of piR-017061 in PDAC, a piRNA that is located within the sno-HBII-296A snoRNA.

Conclusion
This study underlines the role of miRNAs in PDAC and provides evidence for differentially regulated miRNAs that have not been previously implicated in PDAC. Additionally, we provide evidence that novel sncRNA classes, snoRNAs and piRNAs are differentially regulated between normal pancreas and PDAC tissues. Using a bioinformatics approach, we connect mRNA sequencing data with miRNA expression to assign potential functions to miR-802 and other miRNAs. Furthermore, we provide evidence for the differential expression of a variety of lncRNAs in pancreatic cancer.

Methods
Tissue samples were obtained from six patients with PDAC who underwent resection at the Department of Surgery, Technical University of Munich, Germany. Normal pancreatic tissue samples from five patients without pancreatic ductal adenocarcinoma were used as controls.
Tissue collection was approved by the Ethics Committee of the Technical University of Munich and informed consent was obtained from all patients. Tissue were collected directly in the operating theatre and were immediately stored at −80°C until analysis.
Isolation of RNA 20 mg of frozen tissue were disrupted and homogenized (TissueLyser, Qiagen) and RNA was isolated (NucleoSpin miRNA Kit, Macherey-Nagel) in two fractions (small RNA < 200 nt and large RNA > 200 nt).

Cell culture and transfection
The pancreatic cancer cell line MiaPaCa was maintained at 37°C in a humified atmosphere of 5% CO2 and 95% air in Dulbecco's modified Eagle Medium (Life Technologies, Inc, Darmstadt, Germany). The cells were transfected with Lipofectamine2000 (Life Technologies, Inc, Darmstadt, Germany) according to the manufacturer's protocol with either PCMV-MIR-802 or PCMV-MIR-Control (OriGene, Rockville, USA).

RT-qPCR
MiRNA was extracted from MiaPaCa cells with the miR-Neasy Mini Kit (Qiagen, Hilden, Germany). Reverse transcription was performed using the RevertAidH Minus First Strand cDNA Synthesis Kit (Thermo Scientific, Braunschweig, Germany) using specific hsa-miR-802 and as control hsa-miR-16 primer. Amplification of cDNA was performed using the TaqMan Small RNA Assays Applied Biosystems (Life Technologies, Inc, Darmstadt, Germany). The primers for the cDNA synthesis and for the TaqMan analysis were included in the kit from TaqMan Small RNA Assays Applied Biosystems (Life Technologies, Inc, Darmstadt, Germany). For quantification of miRNA-802 expression the results were normalized against miR-16 expression.

Western blot
For quantification of protein levels after miRNA expression, transfected MiaPaCa cells were fractionated using the NE-PER Nuclear and Cytoplasmic Extraction Reagents (Thermo Scientific, Braunschweig, Germany). 20 μg of protein from the nuclear fraction was loaded onto a 10% polyacrylamide gel and was then electrophoretically transferred to a nitrocellulose membrane. The membrane was blocked with Tween-20 (0.05%)-TBS (pH 7.4; 0.1 M Tris Base, 1.4 M NaCl) containing 5% milk, followed by incubation with respective primary antibody α-TCF4 (LS-Bio, Eching, Germany) at a concentration of 1:1000 or as a control α-β-Actin (Abcam, Cambridge, UK) with a concentration of 1:2000. Membranes were washed with Tween-20 (0.05%)-TBS and were incubated with a horseradish peroxidase (HRP)-conjugated anti-rabbit antibody (1:5000). Signals were detected using the enhanced chemiluminescence system (ECL, Amersham Life Science Ltd., Bucks, UK). Films were scanned with a CanoScan 9900F scanner (Canon, Japan). Protein levels were quantified using the Odyssey software LI-COR and normalized to the β-Actin control.

Confirmation of MACE results by qRT-PCR of selected genes
LincRNA and mRNA expression analysis was carried out with the QuantiTect Multiplex PCR Kit (Qiagen) in combination with Superscript III reverse transcriptase (Life Technologies) and PrimeTime qPCR Assays (IDT). For miRNA detection, we used the miRCURY LNA Universal RT microRNA PCR system (Exiqon) according to the recommendations of the manufacturer. Reverse transcription and PCR-amplification for mRNA expression studies were performed with 50 ng of the large total RNA fraction. All quantitative real-time PCR reactions were carried out on the Lightcycler 480 II (Roche). For mRNAs/lincRNAs the expression of housekeeping gene HPRT1 was measured for data normalization, while miR-16 served as endogenous control for miRNA quantification. Differential expression between control and tumor tissues was assessed using the ΔΔC t method, pvalues were calculated with Wilcoxon's rank sum test.

Immunohistochemistry analysis
Immunohistochemistry was performed using the Dako Envision System (Dako Cytomation GmbH, Hamburg, Germany). Consecutive paraffin-embedded tissue sections (3 μm thick) were deparaffinized and rehydrated using routine methods. Antigen retrieval was performed in citrate buffer (pH 6.0; 10 mM citric acid) in a microwave oven for 10 minutes. Endogenous peroxidase activity was quenched by incubation in TBS (pH 7.6; 0.2 M Tris Base; 1.4 M NaCl) containing 3% hydrogen peroxide at room temperature for 10 minutes. After permeabilization with 0.5% TritonX, nonspecific reactivity was blocked with TBS containing 5% BSA. Sections were incubated with the ZEB1 antibody (ZEB1: Atlas Antibodies #AMAb90510 (1:400)) at 4°C overnight followed by incubation with horseradish peroxidase-linked goat anti-mouse antibody, followed by a color-reaction with diaminebenzidine and counterstaining with Mayer's hematoxylin.

Bioinformatics analysis of MACE data
To remove any PCR-bias, all duplicate reads detected by the TrueQuant technology were removed from the raw datasets. The remaining reads were additionally quality trimmed and the poly(A)-tail was clipped off. After raw data processing, reads were aligned to the human genome with novoalign (http://www.novocraft.com). Annotations for genomic mapping positions were derived by the RefSeq annotation track that includes coding genes as well as lncRNAs (http://genome.ucsc.edu/cgi-bin/ hgTables) and only uniquely mapped reads were taken into account. Normalization and test for differential gene expression between normal and tumor tissue were calculated using the DESeq R/Bioconductor package [58]. To account for multiple testing, the false discovery rate (FDR) was estimated. Genes with FDR < 0.05 and |log2fc| > 1.6 were considered as differentially expressed.
Bioinformatics analysis of small RNA-seq data The sRNA-seq data was quantified and tested for differential expression with omiRas [27]. Briefly, for each small RNA-seq library, data processing started with 3′ adapter clipping by a local alignment of the adapter sequence to each read. Subsequently, Illumina's marked quality region was removed and the reads were summarized to UniTags. Singletons were removed from the data set and the remaining tags were mapped to the human genome (hg19) with bowtie. The mapped tags were annotated with the help of various models of coding and non-coding RNAs retrieved from the UCSC table browser. Tags mapping to exonic regions of coding genes were excluded from further analysis. Non-coding RNAs were quantified in each library. For tags mapping to multiple genomic loci the number of reads corresponding to the tag was divided by the number of mapping loci. Normalization and test for differential expression was performed in the same way as described for mRNAs.

Evaluation of normal pancreas function in control libraries
To underline the usability of apparently normal pancreas tissue samples to serve as healthy controls, all genes were sorted in ascending order according to their normalized mean expression in normal sequencing libraries.
To determine genes encoding pancreas specific proteins, functional annotations of the fifty most highly expressed genes were extracted from genecards [59].

Enrichment analysis
Differentially expressed genes were uploaded to the Database for Annotation, Visualization and Integrated Discovery (DAVID) resource [60] using an enrichment cutoff of FDR < 0.05.

Co-expression analysis
A list of all human transcription factors (TFs) was received from AnimalTFDB [61]. Individual expression values of significantly upregulated (in PDAC) transcription factors were clustered using the k-means method with PCC as a distance measure and six initial clusters. The median expression of all transcription factors for each sample was taken as the representative expression for the cluster. To detect the influence of significant miRNA loss on transcription factor upregulation, PCCs were calculated for each TF/TF, miRNA/miRNA and miRNA/TF pair.

Network analysis
Interactions between miRNAs and transcripts with negatively correlated differential expression were detected with omiRas, using a minimum database overlap of three required interaction databases. Additionally, the STRING database v 9.0.5 [62] was used with standard parameters to detect interactions between gene products of differentially expressed genes.