Transcriptional variations in the wider peritumoral tissue environment of pancreatic cancer.

Transcriptional profiling was performed on 452 RNA preparations isolated from various types of pancreatic tissue from tumour patients and healthy donors, with a particular focus on peritumoral samples. Pancreatic ductal adenocarcinomas (PDAC) and cystic tumours were most different in these non-tumorous tissues surrounding them, whereas the actual tumours exhibited rather similar transcript patterns. The environment of cystic tumours was transcriptionally nearly identical to normal pancreas tissue. In contrast, the tissue around PDAC behaved a lot like the tumour, indicating some kind of field defect, while showing far less molecular resemblance to both chronic pancreatitis and healthy tissue. This suggests that the major pathogenic difference between cystic and ductal tumours may be due to their cellular environment rather than the few variations between the tumours. Lack of correlation between DNA methylation and transcript levels makes it unlikely that the observed field defect in the peritumoral tissue of PDAC is controlled to a large extent by such epigenetic regulation. Functionally, a strikingly large number of autophagy-related transcripts was changed in both PDAC and its peritumoral tissue, but not in other pancreatic tumours. A transcription signature of 15 autophagy-related genes was established that permits a prognosis of survival with high accuracy and indicates the role of autophagy in tumour biology.

Pancreatic cancer is one of the most aggressive tumours at all. Despite an incidence that represents only 3% of all cancer cases in industrialised countries, it is the fourth most common cause of tumour‐related deaths in the Western world.1, 2 Pancreatic ductal adenocarcinoma (PDAC) is accounting for >90% of cases and has the worst prognosis; the other tumour types are less lethal. Most PDAC patients die within a year of diagnosis; the overall 5‐year survival rate is about 5%. There is currently no efficient treatment available except surgery, which can only be applied to 10% to 20% of cases, however.3 For an improvement of the clinical situation, prognostic markers are required that allow predicting clinical progression more accurately.4 Several studies on RNA expression in pancreatic cancer have been performed (http://www.pancreasexpression.org/cgi-bin/pancexp/DataSets.pl). Substantial variations have been recorded between different studies, documenting technical variance as well as tumour heterogeneity. Several of these studies dealt with relatively few RNA samples and material from one hospital source only. In addition, only relatively few data are available on pancreatic tumour types other than PDAC. The same is true for non‐tumorous pancreatic tissues and especially so for peritumoral samples. A molecular differentiation of these tissue entities and a detailed understanding of their relation to the actual tumours are missing.

In order to provide a reliable source of information at several molecular levels on a large set of tumour and control samples—including tissue samples that were close to tumour but not part of it—we performed extensive analyses on originally >1,000 pancreatic tissue samples collected by surgery at three major European pancreas clinics. At the DNA level, results on the mutational status of the KRAS and CDKN2A genes and their prognostic significance have been reported, for example.5 Furthermore, suitability of microRNA variations and DNA methylation for diagnosis has been studied6, 7, 8 and followed up by detailed investigations of their functional contributions to the disease9 and the identification of novel routes for therapy.10

Here, we report about an analysis of transcriptional variations at the mRNA level performed on the basis of this large dataset. We identified significant variations between the various tumour forms, but also detected unexpected similarities between mRNA expression patterns. Particularly striking was a substantial degree of similarity in transcriptional regulation between PDAC and the surrounding peritumoral cellular environment, indicating some kind of field defect. Interestingly, the transcriptional variation did not much coincide with changes in the DNA‐methylation levels, which have been implicated in field defects.11, 12 Comparing cystic tumours and PDAC, most transcriptional differences took actually place in this peritumoral environment, while the transcriptional patterns in the tumour tissues were rather similar, suggesting an involvement of the wider cellular environment of a tumour in its pathology. Looking at the data from a functional angle, we inferred relevant pathways and possible functional consequences. Our findings highlight the importance of autophagy‐related transcript expression in the peritumoral environment of pancreatic tumours and the potential role of autophagy‐related genes for prognosis and as legitimate targets for therapeutic intervention schemes.

Material and Methods

Tissue samples and histopathology

Human pancreatic tissue samples were collected during surgery. In all cases, written informed consent was obtained from the patients. The study was approved by the local ethics committees at the universities of Heidelberg, Verona and Liverpool. All “normal” samples were healthy pancreas tissues of donors, who had no pancreatic disease. The samples were snap‐frozen in liquid nitrogen directly after resection and subsequently stored at −80°C. All samples were analysed at DKFZ following identical procedures. The frozen tissue was cut into slices of 15 µm thicknesses with a Leica CM 1850 UV cryotome at −34°C; three slices were picked from the top, middle and bottom third of a tumour and used for histopathology. All remaining slices were mixed to assure equal representation of the entire tissue sample and split into three aliquots, which were used for separate preparations of DNA, RNA and protein.

For a histopathological assessment of each sample's cellular composition, the three tissue sections were stained with hematoxylin and eosin (H&E staining). Selected samples were subjected to Masson trichrome staining (Sigma‐Aldrich, Munich, Germany). They were scanned with a ScanScope GL system (Aperio Technologies, Vista) and visualised using the accompanying software. For each tissue sample, different pathologists evaluated independently the histology and estimated the percentages of normal, tumour and stroma cells.

Immunohistochemistry

Morphological evaluation of the H&E‐stained tissues was supplemented by immunohistochemical staining of the major non‐malignant compartments. Acinar, stromal and immune cells were visualised with antibodies targeting amylase‐A2‐alpha (AMY2A; sheep IgG, abcam #ab18934; Abcam, Camridge, UK), smooth muscle actin alpha (SMACTA2; mIgG2a, DAKO #M 0851; Dako, Hamburg, Germany) and common leucocyte antigen CD45 (mIgG1, DAKO #N1514), respectively. Formalin‐fixed, paraffin‐embedded (FFPE) sections were stained according to a standard protocol.13 In brief, 4 µm‐thick sections were heated to 96°C in citrate buffer (pH 6) for 30 min to retrieve the antigens. They were blocked with methanol containing 3% H2O2 and universal blocking reagent (BioGenex, San Ramon) and exposed to primary antibodies at 4°C overnight. After washing in TBS with 0.05% Tween‐20, slides were exposed for 45 min to anti‐sheep (KPL #5220–0372) or anti‐mouse (DAKO #K400) secondary antibodies labelled with horseradish peroxidase, then incubated with DAB reagent (Dako) for 1 hr and counterstained with hematoxylin. The images were recorded using a light microscope equipped with the AxioVision software (Zeiss, Oberkochen, Germany).

DNA methylation profiling

DNA was isolated with the AllPrep Isolation kit (Qiagen, Hilden, Germany), following the manufacturer's protocol. We performed on 12 normal and 12 PDAC samples the Infinium Human Methylation 450 BeadChip assay of Illumina (San Diego), which interrogates 485,000 methylation sites across the human genome, using 1 µg of DNA per sample. Bisulfite‐converted DNA (EpiTect Bisulfite kit; Qiagen) acted as template for whole‐genome amplification, enzymatic digestion, followed by a DNA clean‐up process and hybridisation to the BeadChip. The samples were washed and scanned with the BeadArray Reader (Illumina). From the signal intensities, the degree of DNA methylation was analysed using the Illumina Bead studio software.

Transcriptional profiling

For RNA isolation, the frozen tissue slices were submerged in liquid nitrogen and gently ground with a polypropylene micropestle (Eppendorf, Hamburg, Germany) in a 2 ml Eppendorf tube. Total RNA was isolated with the AllPrep Isolation kit (Qiagen), following the manufacturer's protocol. RNA integrity was evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto). Only samples with an RNA integrity number of at least seven were used for further analyses. The total RNA prepared from individual samples was analysed on Sentrix Human‐6v3 Whole Genome Expression BeadChips (Sentrix Human WG‐6; Illumina). To synthesize first and second strand cDNA and for amplifying biotinylated cRNA, the Illumina Totalprep RNA Amplification kit was used. Hybridisation to the BeadChip was performed according to the manufacturer's instructions. Subsequently, the arrays were scanned with a BeadArray Reader (Illumina).

Data analysis

Preprocessing and quality control

Raw data were exported from the Illumina Beadstudio software and processed by R/Bioconductor scripts.14 The data was quantile normalised and log2 transformed. Distribution and quality of the expression data was performed using principle component analysis and hierarchical clustering. The raw and normalised data are accessible at the public database ArrayExpress (DNA methylation profiling ID: E‐MTAB‐3855; password “pyzqdbii”; transcriptional profiling ID: E‐MTAB‐1791; password: “rpqqrysi”).

Differential expression analysis

Significant differentially expressed transcript features were detected using the LIMMA package of R/Bioconductor15 by pairwise comparisons of the groups (e.g., PDAC vs. normal; chronic pancreatitis vs. normal; etc.). The resulting p‐values were adjusted for multiple testing using Benjamini‐Hochberg's false discovery rate (FDR) method; features with a FDR < 0.01 and an absolute log2‐fold change |log2FC|>0.5 were considered significant. For a functional enrichment analysis, the Ingenuity PA software tool (Ingenuity Systems, http://www.ingenuity.com) was applied.

Survival analysis

Cox's proportional hazards models were used for exploring the relationship between the survival of a patient and several explanatory variables, including phenotypic parameters and gene expression. Analysis was performed on PDAC patients utilising the survival package of R/Bioconductor. This analysis returns statistical significance of prognostic variables included into the model and allowed to estimate the risk of death for individuals. A positive regression coefficient (or log hazard ratio) for an explanatory variable at high significance level means that the hazard is higher, and thus the prognosis is worse. Conversely, a negative regression coefficient implies a better prognosis for patients. We started with checking the effect of phenotypical parameters: age, gender, cancer stage, smoking and alcohol intake. Then the effect of gene expression was investigated for each gene independently and FDR‐adjusted p‐values of a Wald test were assigned to them.

Signature selection amongst autophagy‐related genes

We used the intersection among significant genes in order to identify autophagy‐related genes, which can be used as prognostic signature for PDAC. Before analysis, we summarized microarray features to the gene level in order to merge data from features targeting the same gene. Genes were selected that were differentially expressed, linked to survival by Cox regression (FDR < 0.01) and annotated as autophagy‐related genes in relevant literature and databases16, 17 (see also Supporting Information Table S1). In order to transform expression intensity values of genes to a single value, we used an approach presented before.18 Gene expression values were median‐centred, and expression values of genes with a negative Cox's coefficient were inverted by multiplication by minus one in order to account for a survival effect. Then, gene expression of each patient was summed up forming a single score. Predictive power of the signature was characterised using support vector machine (SVM) as a classifier. Cancer patients were divided into two groups with “good” and “poor” prognosis based on their survival in relation to median survival time of the patients.19 Gene expression was used as input to a support vector machine classification; output was considered as a binary signal for “good” or “poor” survival. During performance testing, 80% of the patients were randomly selected as a training set, keeping a correct proportion of good/poor prognosis patients. This operation was repeated 1,000 times in order to characterise confidence intervals for accuracy (ratio between numbers of correctly classified patients to total number).

Pathway analysis of differentially expressed genes

Differentially expressed transcripts were mapped onto a molecular network developed from information contained in the Ingenuity knowledge base (Ingenuity Systems). Gene networks were generated based on their interconnectivity. Ingenuity Pathway Analysis ranks the resulting networks by calculating a significance score corresponding to the negative log of p values. Furthermore, pathway core analysis identified the pre‐specified canonical pathways that were most over‐represented in the data set. Fisher's exact test was used to calculate a p values for the association between the genes in the data set and the canonical pathway or network.

RT‐PCR confirmation

For RT‐PCR, 1 μg total RNA of 20 samples per group (normal pancreas, PDAC and PDAC macro‐environment) was reverse‐transcribed using the ProtoScript M‐MuLV First Strand cDNA Kit (New England Biolabs, Frankfurt, Germany). Quantitative RT‐PCR was performed in triplicate on a LightCycler 480 (Roche Diagnostics, Mannheim, Germany) with a pre‐amplification incubation of 95°C for 10 min, followed by 45 cycles of 95°C for 10 sec and 55°C for 30 sec. The following molecules of the QuantiTect Primer Assays (Qiagen) were used as primers: Hs_ACTB_1_SG, Hs_PKM_1_SG, Hs_PLK4_1_SG, Hs_PRAF2_1_SG and Hs_BUD13_1_SG. BUD13 was used as reference gene as it showed constant mRNA expression in all 452 analysed pancreatic tissue samples and an adequate expression level. Data were analysed using the LightCycler software (Roche).

Results

Comparison of transcriptional variations in different tissue types

The study was performed on initially 1,004 pancreatic tissue samples collected at three clinics in Heidelberg (Germany), Liverpool (UK) and Verona (Italy) by resection from cancer patients or donors who had pancreatic tissue removed for reasons other than cancer. Molecular sample analysis was performed in one central laboratory, following the same protocols throughout. From each sample, three tissue slices were evaluated by experienced pathologists in a central pathology review prior to further analysis, estimating the percentages of normal, tumour and stromal cells as well as the degree of inflammatory infiltration. In total, 452 RNA preparations met our quality standards (see Methods section) and were included in the transcriptional profiling analysis. The quality assessment was done prior to any subsequent analysis, thus purely based on RNA quality and unbiased. The percentage of only about 45% good‐quality preparations is not surprising when dealing with clinical material, and particularly with pancreas samples, which contain a very high level of RNases. Many samples of lower quality could have been used on the microarrays. However, the significantly higher background combined with weaker signal intensities that results from degraded RNA would have made the data far less accurate and reproducible and could not have been offset by the then larger number of samples in the overall analysis.

The distribution of the analysed samples was studied by a principle component analysis (Fig. 1), already indicating differences between the tissue types. In addition, few apparent discrepancies of the histological and molecular stratification could be observed. Of particular interest were samples of non‐tumorous tissue that had been located next to the actual tumour. These peritumoral tissues had been in a distance of about 1 to 10 mm from the tumour and did not exhibit any tumour cell content in the histochemical analyses. They consisted of parenchyma and stroma and are referred to as “macro‐environment” below. The 452 high‐quality samples represented 195 cases of PDAC, 30 cases of PDAC macro‐environment (N_PDAC), 24 cystic tumours (CT), 22 macro‐environmental tissues from next to cystic tumours (N_CT), 59 samples of chronic pancreatitis (CP), and 41 healthy pancreatic tissues from non‐cancer patients (N). The subtypes of the cystic tumours and the relevant macro‐environments are listed in Supporting Information Table S2. Other neoplasms, for which RNA was isolated but ignored in the analysis reported below, were 18 endocrine tumours, 2 macro‐environmental samples of endocrine tumours, 31 other pancreatic tumours, 15 related macro‐environment tissues as well as 15 tissues from the macro‐environment of CP (N_CP). Information about clinical patient parameters is given in Table 1.

Figure 1

Principle component analysis of the samples based on their transcript profiles. The colour code of the tumour types is given at the top. The plot shows a high degree of similarity of PDAC and cystic tumours, indicates a distinct difference between their macro‐environments and highlights a similarity of the macro‐environments of cystic tumours and healthy tissues.

Table 1

Clinical characteristics of patient cohort

	N	CP	PDAC	N_PDAC	CT	N_CT	Others
No. of patients	41	58	195	30	24	22	82
Gender (male/female)	26/15	48/10	109/86	21/9	7/16	4/18	50/32
Age at surgery, median (range)	46.0 (16–74)	47.1 (13–73)	63.4 (40–85)	60.8 (34–84)	62.0 (23–75)	57.1 (38–75)	55.7 (13–86)
Stage
0	n/a	n/a	–	–		–	–
IA	n/a	n/a	–	–	2	–	1
IB	n/a	n/a	1	–	0	–	1
IIA	n/a	n/a	21	5	3	3	5
IIB	n/a	n/a	123	18	6	1	24
III	n/a	n/a	7	–	–	–	2
IV	n/a	n/a	17	2	1	–	4
Median survival‐time in months, (range)	n/a	n/a	24.7 (1–159)	19.7 (1–65)	15.9 (1–36)	48.5 (1–141)	18.22 (1–54)

Details are listed of the clinical parameters of the patients from whom the 452 RNA‐preparations were isolated and subsequently analysed. N: healthy tissue; CP: chronic pancreatitis; PDAC: pancreatic ductal adenocarcinoma; N_PDAC: macro‐environment of PDAC; CT: cystic tumour; N_CT: macro‐environment of cystic tumour; n/a: not applicable

Many genes exhibited changes at the transcript level in the various tissue types compared to samples from healthy donors (Supporting Information Fig. S1; Supporting Information Table S3). Also, a large number of changes were found that were common between different tissues (Fig. 2). It is noteworthy that this means a variation in the same direction (up or down, respectively) as compared to the normal tissue and not an increase in one and a decrease in the other tissue. As a matter of fact, an analysis revealed that basically all genes, which showed a significant variation in two tissues, were similarly regulated either up or down in both when comparing N_PDAC, PDAC and cystic tumours (Fig. 2 e). This lack of inverse regulation suggests that functionally similar cellular effects are triggered by these common transcriptional changes.

Figure 2

Tissue specificity of mRNA level variations. For each tissue type, the number of mRNAs is shown that were significantly differentially expressed in comparison to normal pancreas tissue (N). The numbers in overlap regions stand for genes, regulated similarly in the relevant tissues. (a) Results are presented for PDAC, the related macro‐environment (N_PDAC) and chronic pancreatitis (CP), marked in red, green and yellow, respectively. (b) The panel presents the same for cystic tumours (TC; brown), the related macro‐environment (N_CT; blue) and again chronic pancreatitis (CP; yellow). (c) The macro‐environment of cystic tumours (N_CT) exhibited relatively few variations at the mRNA level that were specific. (d) Presentation of the result of a comparison of all five data sets. (e) Correlation in the direction of variation observed for N_PDAC versus N (top panel) or CT versus N (bottom panel), respectively, in comparison to PDAC versus N. Both axes represent the score shown above the panels, thus focussing on the most significant variations (shown in blue). Grey dots, mostly close to the centroid, represent insignificant changes. All actual data is accessible in Supporting Information Table S3.

Particularly interesting results were obtained from the macro‐environments of PDAC and cystic tumours, which histologically showed a similar cell composition (Supporting Information Figs. S2 and S3). To confirm this further, we also performed immunohistochemical analyses. Acinar, stromal and immune cells were visualised in formalin‐fixed, paraffin‐embedded sections with antibodies targeting amylase‐A2‐alpha, smooth muscle actin alpha and CD45, respectively (Supporting Information Fig. S4). Non‐tumorous peritumoral samples exhibited no abnormalities, such as acinar loss, tumour cell invasion, stromal hyperplasia or inflammatory infiltrates. Only two of the PDAC macro‐environments exhibited a GGT > GAT mutation of KRAS in codon 12; all other samples had no KRAS mutations. Compared to the transcript profiles obtained from normal tissue, a large number of genes were found differentially transcribed in PDAC macro‐environment, indicating that the phenotypically non‐tumorous appearance does not represent the molecular status (Fig. 2 a). In addition, there was a substantial overlap of 2,641 genes with the results obtained from PDAC, of which 1,997 were also changed in CP. At a transcriptional level, the PDAC macro‐environment and tissue of patients with CP were similarly different to normal tissue. However, PDAC macro‐environment exhibited significantly more molecular resemblance with tumour tissue than with CP. For cystic tumours, there was a rather different picture. The overlap with the transcriptional pattern of CP was identical in number to the overlap of PDAC and CP (Fig. 2 b). However, in contrast to the PDAC result, the macro‐environment of cystic tumours behaved very similar to normal pancreatic tissue with only 343 differentially transcribed genes as opposed to 2,909 genes in PDAC macro‐environment. In a comparison of the macro‐environments of PDAC and cystic tumours (Fig. 2 c), only 45 differentially expressed transcripts were found to be specific for the macro‐environment of cystic tumours. All this is also graphically represented in the principle component analysis (Fig. 1).

Variations that are specific to PDAC and its macro‐environment

Since PDAC represents the vast majority of clinical pancreatic cancer cases and has the worst prognosis, we focussed our analysis on this tumour type. More than half of the 5,196 genes that exhibited significant variations in PDAC compared to healthy tissue were also differentially expressed in PDAC macro‐environment or CP tissues (Fig. 2 a). The 1,997 regulated genes shared between the three tissue types mostly represent changes that are associated with inflammation. The most overrepresented canonical pathways defined by these genes are relevant for immunological and inflammatory response (Supporting Information Fig. S5). Next to the shared transcript variations, there were (marker) genes that exhibited expression changes that were unique to each tissue type. The six most significant biological functions associated with the 2,373 unique PDAC expression markers are: cellular growth and proliferation; cellular movement; cell death and survival; cancer; cell cycle; as well as organismal injury and abnormalities (Fig. 3). There were no over‐represented functions that are associated with inflammation. In contrast to this, the 127 and 162 markers that were unique to the PDAC macro‐environment or CP, respectively, are genes over‐representing functions associated with both inflammation and cancer. An analysis of the CT expression markers produced results rather similar to that of PDAC (Fig. 3).

Figure 3

Most overrepresented biological functions associated with the 2,373 and 2,619 unique marker genes of PDAC and CT, respectively. Colour‐coded maps of functional predictions resulting from an Ingenuity Pathway Analysis are shown. Each rectangle represents one function. The intensity of the purple colour of a square is proportional to the number of genes that are associated with the respective function. The size of a square reflects the associated negative log10 of the assigned p values. Larger squares indicate a more significant overlap between the genes perturbed in the dataset and the respective function. The top 10 biological functions were ordered in a bar blot according to their significance by the negative log10 of the assigned p values. The complete lists can be found at Supporting Information Fig. 6.

Of particular interest are the 644 regulated transcripts shared between PDAC and its macro‐environment but not regulated in CP (Fig. 2 a; Supporting Information Table S4). The fact that the non‐tumorous macro‐environment tissue exhibited a large number of variations, which are common with PDAC only, suggests the possibility of a cancer field defect at the transcript level. A cancer field defect is defined as a biological cancerisation process in which tissue in relatively large areas beyond the actual tumour is exhibiting epigenetic changes similar to the ones of the actual tumour.11, 12 Methylation changes are expected to result in variations at the transcript level, too. To reveal the degree by which promoter methylation may be responsible for the observed variations at the RNA level, we analysed the genomic DNA of 24 randomly selected samples, 12 each from healthy donors and PDAC patients, from which we had also obtained mRNA profiles. In the set of 644 transcripts shared by PDAC and its macro‐environment, 154 genes exhibited an inverse correlation of promoter methylation and mRNA expression: 115 hyper‐methylated promoters could be linked to down‐regulation of the mRNA level, 39 hypo‐methylated promoters coincided with transcripts that were present in higher abundance. However, for the majority of genes, 490, there was no such correlation of promoter DNA methylation and gene expression. For 76 genes, there was even a concurrent increase or decrease, respectively, of both methylation and mRNA level.

Pathways and gene networks affected in the macro‐environment of PDAC

Transcripts showing significant expression differences in the macro‐environment of PDAC as compared to healthy control samples were submitted to a functional bioinformatics examination using Ingenuity pathways analysis. In this analysis, one particular network of genes was identified. It consists of 73 directly linked genes (Supporting Information Fig. S7; Supporting Information Table S5). Many genes show a gradual change of expression from healthy via CP to PDAC macro‐environment and finally PDAC. Two typical examples are the chemokine receptor CXCR4 (chemokine (C‐X‐C motif) receptor 4) and NUPR1 (nuclear protein transcriptional regulator 1). CXCR4 is a prognostic marker in various types of cancer and a biomarker of migrating pancreatic cancer‐initiating cells in mice.20 It was up‐regulated compared to normal tissue by a factor of 3.51 in CP, 4.38 in the macro‐environment and 6.68 in PDAC. NUPR1 interacts with numerous partners to regulate cell cycle, apoptosis, autophagy, chromatin accessibility, and transcription. Reduced expression promotes pancreatic cancer development.21 It was downregulated by −1.80 (CP), −2.10 (macro‐environment) and −2.41 (PDAC). Numerous genes within the network are essential in Pancreatic Adenocarcinoma Signalling, as well as the Integrin‐Linked Kinase (ILK) and mTOR Signalling pathways. ILK is known to transmit mechanical stimuli to the mammalian target of the mTOR Signalling pathway. mTOR is an evolutionarily conserved protein kinase that integrates both intracellular and extracellular signals and serves as a central regulator of cell metabolism, growth, proliferation, survival, and autophagy. The Akt/mTOR pathway mediates oncogenesis and controls tumour cell growth.22

Based on these findings, we took a closer look at these pathways (Supporting Information Fig. S8). Interestingly, in the Pancreatic Adenocarcinoma Signalling pathway, 35 of the 106 genes (33%) are annotated as relevant to autophagy. For ILK/mTOR Signalling, the absolute number of autophagy‐related genes was even higher with 54 molecules out of 186, although the percentage was slightly lower with 29%. The considerable overrepresentation of autophagy‐associated genes suggested a strong link between autophagy and pancreatic cancer. This was corroborated by the overall number of autophagy‐related genes found to be regulated in PDAC. Of the genes that were assayed in the transcriptional profiling experiments, 512 are annotated to be associated with autophagy processes according to two public autophagy databases.16, 17 Of these, 88 genes were regulated in CP compared to normal tissue. In the macro‐environment, this number increased to 108 genes, while 208 genes were regulated significantly in PDAC (Supporting Information Table S6). Given the obvious importance of autophagy for pancreatic cancer and the fact that about 40% of all autophagy‐related genes were significantly regulated in PDAC or in the PDAC macro‐environment, we focussed further analyses on those genes.

Autophagy‐related genes in PDAC, cystic tumours and macro‐environment

In our data set, there were several autophagy‐related genes that are known markers for pancreatic cancer. Examples are ATG3, ATG4A, ATG4B, ATG4C and ATG5 23 showing down‐regulation in PDAC and its macro‐environment. Also the ductal transcription factor HNF6 (ONECUT1) was less expressed. Such a reduction in HNF6 expression correlates with human pancreatic cancer progression.24 HNF6 has also been described as biomarker for acinar‐to‐ductal metaplasia25 suggesting that a phenotypic switch converting pancreatic acinar cells to duct‐like cells could lead to pancreatic intraepithelial neoplasia and eventually to invasive PDAC. Furthermore, HNF6 is a transcription regulator of the UDP glucuronosyltransferase (UGT) family members. Interestingly, variations in the UGT genotype are associated with an altered risk to pancreatic cancer.26 None of the above genes was regulated in the macro‐environment of cystic tumours. As mentioned before, many mRNAs were significantly differentially transcribed in CP, PDAC macro‐environment, PDAC and cystic tumours, but not in the macro‐environment of cystic tumours (Fig. 2 c). As cystic tumours are known to have a much better prognosis than PDAC,27 factors that are regulated in the macro‐environment of PDAC but not in the macro‐environment of cystic tumours might be relevant for the much more aggressive nature of PDAC.

Autophagy‐related prognostic markers

Given the apparent importance of autophagy‐related genes as specific markers of PDAC and its macro‐environment, we wondered whether they could act as a clinically relevant surrogate for disease and have prognostic significance. To investigate this, the relationship was explored between patient survival and several explanatory variables. First, a set of univariable survival analyses was performed to address the potential influence of the factors age, gender, tumour stage, treatment, smoking, and alcohol consumption. Only tumour stage showed a significant linkage to survival (Wald p values of 0.0013). More specifically, we observed the difference between two groups, early stage (S2) and late stage (S3, S4) patients, with a Cox regression coefficient or log hazard ratio of 1.02 ±0.51 (95% confidence interval). Then, the effect of gene expression was investigated for each gene individually. In total, 35 autophagy‐related transcripts were differentially transcribed in PDAC and statistically linked to survival time (Supporting Information Table S7). Eighteen genes were regulated in both PDAC and its macro‐environment (Table 2). The top‐candidates were PRAF2, PLK4, ACTB and PKM2. The transcript levels were confirmed by qRT‐PCR, fitting to the microarray data in all cases (Supporting Information Fig. S9). Interestingly, neither PRAF2 nor PLK4 had been connected to PDAC so far. With respect to prognosis, the expression of the two genes corresponded well with survival time and matched the performance of hypoxia‐inducible factor 1 alpha (HIF1A; Fig. 4), whose change in expression has been shown to have a strong impact on the prognosis of patients with PDAC.28

Table 2

Result of a Cox regression of 18 autophagy transcripts statistically linked to survival time

Gene	Cox Coeff.	Cox FDR	PDAC log2FC	PDAC FDR	N_PDAC log2FC	N_PDAC FDR
ACTB	1.087	0.0023	1.298	0.000	0.768	0.000
ANTXR1	0.601	0.0088	1.373	0.000	0.843	0.000
CAMK1G	−0.410	0.0008	1.140	0.000	1.015	0.001
DLG4	0.570	0.0001	0.761	0.000	0.902	0.000
DNAJB9	−0.598	0.0069	−1.300	0.000	−0.904	0.000
EIF2AK3	0.846	0.0014	−0.715	0.000	−0.501	0.001
ITPR1	−0.704	0.0017	0.741	0.000	0.816	0.000
MAP2K7	0.697	0.0179	−0.901	0.000	−0.608	0.000
MPDZ	−0.685	0.0038	0.613	0.000	0.729	0.000
MYO5C	−0.407	0.0036	−0.818	0.000	−0.780	0.001
NLE1	1.153	0.0050	−0.728	0.000	−0.532	0.000
P4HB	0.550	0.0057	−1.079	0.000	−0.592	0.003
PKM2	0.617	0.0004	1.475	0.000	0.702	0.000
PLK4	−0.633	0.0034	0.991	0.000	0.529	0.000
PRAF2	1.068	0.0000	0.588	0.000	0.595	0.000
SH3GLB2	−0.550	0.0198	−0.726	0.000	−0.713	0.000
VIM	0.421	0.0166	0.846	0.000	0.865	0.000
WDFY2	−0.822	0.0052	−1.192	0.000	−0.915	0.000

A positive coefficient indicates a worse prognosis, a negative coefficient a protective effect. In addition, the regulation in PDAC and PDAC macro‐environment (N_PDAC) as compared to healthy tissue is shown. FRD: false discovery rate; log2FC: logarithm of fold change

Figure 4

Linkage of gene expression levels in PDAC and patient survival time. Four typical Kaplan–Meyer curves with 95% confidence (dotted lines) are shown. An increase in expression of PRAF2 is linked to poor survival, whereas stronger PLK4 expression predicts better survival. In the third panel, the result is shown for the established prognostic marker HIF1A, which is linked to poor survival. The expression of ABL1 is not correlated to survival time at all and shown as a reference. The p‐values shown are based on a Cox regression of continuous log2 gene expression.

Based on the prognostic value of individual genes and by cross‐referencing this with the list of autophagy‐related genes that are differentially expressed in PDAC and its macro‐environment, we identified the most‐informative signature, which consists of 15 genes: ACTB, ANTXR1, CAMK1G, DLG4, DNAJB9, EIF2AK3, ITPR1, MPDZ, MYO5C, NLE1, P4HB, PKM2, PLK4, PRAF2 and WDFY2 (Supporting Information Fig. S10). Expression of these genes was combined into a single score as described in the Methods section. As a result, the score showed a strong statistical linkage to survival in a Cox regression model (Wald p values = 7 × 10−14, Cox coefficient = 0.11 ±0.03). We observed a statistical link between the calculated score and tumour stage (ANOVA p values of 0.013). Next, we classified patients by a support vector machine as described in the Methods section: our classifier discriminated “good” and “poor” prognosis patients. As a basis, they were divided into two groups based on their survival in relation to the median survival time. The classification was performed 1,000 times using cross‐validation approach and resulted in an average accuracy of 82.6% ±0.1% on training data and 65.4% ±0.5% on test data. The discriminating power of the signature is shown in a Kaplan–Meier plot (Fig. 5), indicating the substantial difference in survival between the patients with low and high score calculated on that signature.

Figure 5

Prognosis of patient survival. Kaplan–Meier curves were calculated based on an expression signature in PDAC of the 15 genes named in the figure (Wald p values of Cox model = 7 × 10−14). In blue, the survival of the patients with low score (below median) is shown; the red line represents the result of the patients with high score (above median). The dotted lines correspond to 95% confidence intervals.

Discussion

We performed transcriptional profiling on various types of pancreatic tissue samples, focussing particularly on changes in tissues that were located next to actual tumours. Comparison of the expression data documented that PDAC and cystic tumours are not that different at the transcriptional level, even though 892 and 516 genes, respectively, had expression patterns that were unique to either tumour type. However, variation in 4,179 genes was in common. In addition, both tumour types shared a rather similar set of transcriptional variations with CP tissues, indicating the substantial contribution of inflammatory aspects that are relevant for tumour pathology. The most significant difference between PDAC and cystic tumours was actually found in their surrounding, non‐tumorous tissues. Histologically, they exhibited a similar cell composition. However, while the PDAC macro‐environment behaved in part like the actual tumour, the macro‐environment of cystic tumours was transcriptionally nearly identical to normal pancreas tissue. This suggests that major differences between cystic tumours and PDAC may not be solely intrinsic to the actual tumours, but could be triggered indirectly by the way tumours influence or are influenced by their wider cellular environment.

Assuming that a field defect—an influence of the tumour on the adjacent tissue11, 12—could be responsible for the large number of genes, which were identically regulated in PDAC and its macro‐environment, and that it may be resulting from DNA methylation processes, which have been implicated in field defects, we studied how the degree of promoter methylation correlated with the transcriptional variations. Surprisingly, the degree of inverse correlation—hyper‐methylation and low expression, hypo‐methylation and high expression—was rather low overall. One possible explanation for this observation could be the presence of processes other than DNA methylation, which regulate transcription. Any such mechanism must be rather efficient in transporting the information, since the effect could be detected over distances of several millimetres. MicroRNA transport in cellular vesicles, such as exosomes, and absorption by recipient cells has been shown to have such effects.29 In the intracellular space, one would expect a gradient in the concentration of exosomes radiated by tumours and a related gradient in the degree of cellular transcriptional variations. The actual distance of the macro‐environmental samples analysed in this study was inadequately annotated so as to identify such an effect.

Alternatively, it could be that only very few genes and their methylation may be initiating the field defect. Gene NR5A2 (encoding the nuclear receptor subfamily 5, group A, member 2) was found to be one of the genes, which was both hyper‐methylated and down‐regulated in PDAC and PDAC macro‐environment. NR5A2 had not been described as a methylation marker of PDAC before. However, single nucleotide polymorphisms in the vicinity of NR5A2 have been linked through genome‐wide association studies to the risk of developing PDAC, suggesting a broader role of this gene in pancreatic homeostasis and disease.30 It was shown that NR5A2 heterozygosity correlates with pancreatic damage in the progression of mutant KRAS‐driven preneoplastic lesions, suggesting that NR5A2 could contribute to PDAC through its role in the recovery from pancreatitis‐induced damage.31 The gene might be an interesting candidate for exploring a potential cancer field defect in pancreatic cancer. With respect to clinical utility, a field defect in the non‐tumorous tissue and the identification and validation of relevant molecular variations at the RNA level or the degree of DNA methylation could allow the establishment of effective disease biopsy markers, since they would also be present in a wider distance to the actual tumour and therefore be easier to collect.

One cannot rule out that other mechanisms than the ones above might be responsible for the field defect. For example, surgeons and pathologists observe that peritumoral tissues in patients with PDAC are apparently altered phenotypically, while tissues adjacent to cystic or neuroendocrine tumours often look quite normal. The mere degree of expansion of PDAC and the resulting obstruction of the macro‐environment could be triggering transcriptional variations that also occur in PDAC tissue but not in or around the slowly growing cystic tumours, thus providing a simple mechanical explanation for the differences. Further studies are required to define the molecular processes in more detail.

Our findings particularly indicate the relevance of autophagy‐related transcripts in the macro‐environment of pancreatic tumours and their potential role as prognostic markers. Autophagy is an evolutionarily conserved catabolic process, by which a cell digests its own cytoplasmic content. Autophagy is activated in reaction to multiple stress factors during cancer progression, such as hypoxia and poor nutrient supply32 and has an important role in tumour development.33 Overall, the dynamic role of autophagy in cancer appears to be complex and context‐dependent. On the one hand, it could function as a tumour suppressor, whose inactivation promotes tumorigenesis,34 on the other hand, it may act as a pro‐survival pathway that helps tumour cells to handle metabolic stress and to resist chemotherapeutic agents.35 In pancreatic cancer, autophagy is required for tumour growth36 and mediates survival of pancreatic‐tumour‐initiating cells in a hypoxic microenvironment.37 Inhibition of autophagy reduced pancreatic cancer growth independent of the p53 status.36 Furthermore, autophagy was shown to be essential for oncogenic KRAS‐induced malignant cell transformation. There are several early‐phase clinical trials in progress targeting the autophagic machinery, among them a study of MEK1/2 and AKT inhibitors in patients with KRAS‐driven pancreatic tumours.38

Prognostic significance of autophagy‐related protein expression in resected pancreatic ductal adenocarcinoma has been described.39 The study reports about a correlation of the expression of five proteins with disease prognosis and was based on immunohistochemistry on 73 tumour tissue sections. Our analysis studied all genes and has a focus on the variations that occur in the macro‐environment of the tumour. The regulation of several autophagy‐associated genes correlated with survival‐time. Of special interest are the prognostic factors that were particularly regulated in pancreatic cancer and its macro‐environment, such as PRAF2 and PLK4, as these transcripts might be useful in distinguishing inflamed versus oncogenically transformed regions of the pancreas. PRAF2 (Prenylated Rab acceptor 1 domain family, member 2) is a small transmembrane protein with a putative role in transport from the endoplasmic reticulum to the Golgi apparatus. It induces apoptotic cell death upon expression and is counteracted by Bcl‐xL.40 It stimulates cell proliferation and migration and predicts poor prognosis in neuroblastoma and glioma.41, 42 In a genome‐wide siRNA screen aiming at the identification of new genes involved in autophagy, PRAF2 was found as a gene mediating autophagy.43 PLK4 (Polo‐like kinase 4) is a conserved upstream regulator of centriole duplication. It is aberrantly expressed in different tumour types,44 causing a loss of centrosome numeral integrity, thereby promoting genomic instability. PLK4 interacts with Cep63, which was shown to be degraded in autophagy‐controlled centrosome number variations.45 Recently, a cell line model of centrosome amplification was described, which exhibited elevated levels of autophagy upon induced expression of PLK4.46 A drug discovery programme identified a potent and selective small molecule inhibitor of PLK4,47 which may have therapeutic implications for PDAC.

Already individual differentially expressed autophagy‐related transcripts allowed a disease prognosis that is equivalent to markers reported before. By combining 15 autophagy‐related genes, we obtained a signature, which allowed an even better classification. Even though some transcripts exhibited an only modestly differential expression in PDAC, they documented a significant prognostic power. These transcripts may offer a means for a better prognosis after tumour resection, although limitations apply for utilising actuarial probabilities for a prediction of the actual survival of an individual patient.48 However, the analysis highlights the importance of autophagy for tumour pathology and indicates that this process is likely to be highly relevant for future treatment strategies and monitoring of their success.

Author Contributions

A.S.B., A.S., M.W.B., J.P.N. and J.D.H. conceived the study; A.S.B., N.G., S.B., A.H. and M.B. performed experiments; A.S.B., P.V.N., N.G., W.G., E.C., A.M., O.S., T.H., L.V., A.S. and S.K. analysed and interpreted data; N.G., O.S., T.H., A.S., M.W.B. and J.P.N. provided vital reagents; A.S.B., P.V.N., N.G. and J.D.H. wrote the manuscript; all authors contributed to manuscript writing and editing.

Supporting Information 1

Click here for additional data file.

Supporting Information 2

Click here for additional data file.

Supporting Information 3

Click here for additional data file.

Supporting Information 4

Click here for additional data file.

Supporting Information 5

Click here for additional data file.

Supporting Information 6

Click here for additional data file.

Supporting Information 7

Click here for additional data file.

Supporting Information 8

Click here for additional data file.

Supporting Information 9

Click here for additional data file.

Supporting Information 10

Click here for additional data file.

Supporting Information 11

Click here for additional data file.

Supporting Information 12

Click here for additional data file.

Supporting Information 13

Click here for additional data file.

Supporting Information 14

Click here for additional data file.

Supporting Information 15

Click here for additional data file.

Supporting Information 16

Click here for additional data file.

Supporting Information 17

Click here for additional data file.

Supporting Information 18

Click here for additional data file.