Differences between Well-Differentiated Neuroendocrine Tumors and Ductal Adenocarcinomas of the Pancreas Assessed by Multi-Omics Profiling

Most pancreatic neuroendocrine tumors (PNETs) are indolent, while pancreatic ductal adenocarcinomas (PDACs) are particularly aggressive. To elucidate the basis for this difference and to establish the biomarkers, by using the deep sequencing, we analyzed somatic variants across coding regions of 409 cancer genes and measured mRNA/miRNA expression in nine PNETs, eight PDACs, and four intestinal neuroendocrine tumors (INETs). There were 153 unique somatic variants considered pathogenic or likely pathogenic, found in 50, 57, and 24 genes in PDACs, PNETs, and INETs, respectively. Ten and 11 genes contained a pathogenic mutation in at least one sample of all tumor types and in PDACs and PNETs, respectively, while 28, 34, and 11 genes were found to be mutated exclusively in PDACs, PNETs, and INETs, respectively. The mRNA and miRNA transcriptomes of PDACs and NETs were distinct: from 54 to 1659 differentially expressed mRNAs and from 117 to 250 differentially expressed miRNAs exhibited high discrimination ability and resulted in models with an area under the receiver operating characteristics curve (AUC-ROC) >0.9 for both miRNA and mRNA. Given the miRNAs high stability, we proposed exploring that class of RNA as new pancreatic tumor biomarkers.


Introduction
Pancreatic neuroendocrine tumors (PNETs) arise from the endocrine portion of the pancreas, and pancreatic ductal adenocarcinomas (PDACs) originate from epithelial cells. Although most are indolent, PNETs are unpredictable, ranging from nearly benign to aggressive, metastasizing neoplasms.

Mutational Status
To compare potential driver mutations in neuroendocrine tumors and ductal adenocarcinomas of the pancreas, targeted sequencing of 409 genes was performed in 13 NETs (nine PNETs and four INETs) and 8 PDACs; detailed characteristics of the patients are shown in Table 1. The median mapped read count was 20,820,942. The median percentage of bases with greater than 100× coverage was 98.41%. Each read was assigned to its respective amplicon to perform mutation calling and estimate the mutant allele frequency.
After quality filtering, we identified 1288 unique variants across 234 genes, ranging from 106 to 206 variants for each tumor sample. Among them, 153 were considered pathogenic or likely pathogenic according to the VAST or CHASM algorithm or had a high impact on protein coding (stop-gained, frameshift, etc.) (Table S1). Of these, variants likely to be deleterious were found in 50, 57, and 24 genes in PDACs, PNETs, and INETs, respectively. Ten genes and 11 genes contained a pathogenic or likely pathogenic mutation in at least in one sample of all tumor types and in PDACs and PNETs, respectively, while 28 genes, 34 genes, and 11 genes were found to be mutated only in PDACs, PNETs, and INETs, respectively ( Figure 1, listed in Table 2).  The seven most frequently mutated genes (HNF1A, CREBBP, RNF213, KRAS, CASC5, APC, and TP53) were found in 50% or more of PDAC samples. Two genes (CASC5, APC) and six genes (HNF1A, CREBBP, KAT6B, MLL2, MLL3, and JAK3) were found to be mutated in 44% and 33% of PNET samples, respectively. Eleven genes (STK11, CASC5, CREBBP, RNF213, PAX7, NFKB2, FLT4, NOTCH1, FLT3, ZNF521, and EP300) were found to be mutated in two of four INET samples. Of these, the most frequently mutated genes were CREBBP and CASC5 in all three tumor types, HNF1A  The seven most frequently mutated genes (HNF1A, CREBBP, RNF213, KRAS, CASC5, APC, and TP53) were found in 50% or more of PDAC samples. Two genes (CASC5, APC) and six genes (HNF1A, CREBBP, KAT6B, MLL2, MLL3, and JAK3) were found to be mutated in 44% and 33% of PNET samples, respectively. Eleven genes (STK11, CASC5, CREBBP, RNF213, PAX7, NFKB2, FLT4, NOTCH1, FLT3, ZNF521, and EP300) were found to be mutated in two of four INET samples. Of these, the most frequently mutated genes were CREBBP and CASC5 in all three tumor types, HNF1A and APC in both PDACs and PNETs, while two, four, and eight genes were mutated only in PDACs, PNETS, and INETs, respectively ( Figure 2, Table 3). and APC in both PDACs and PNETs, while two, four, and eight genes were mutated only in PDACs, PNETS, and INETs, respectively ( Figure 2, Table 3).
Of genes previously reported as commonly mutated in PNETs, only PDGFRA and MTOR were found to be mutated in a few tumor samples, and other genes, including ATM, ATRX, TSC2, DAXX, VHL, and MEN1, were mutated only in individual samples.

Gene Expression Analyses
We analyzed gene expression in tissue samples by deep sequencing of the mRNA and miRNA transcriptomes. Differentially expressed genes were identified using Wald's test, and principal component analysis (PCA) of these genes revealed visible separation of PDACS from NETs and some, but not so clear, second principal component separation of PNETs from INETs ( Figure 3).  Of genes previously reported as commonly mutated in PNETs, only PDGFRA and MTOR were found to be mutated in a few tumor samples, and other genes, including ATM, ATRX, TSC2, DAXX, VHL, and MEN1, were mutated only in individual samples.

Gene Expression Analyses
We analyzed gene expression in tissue samples by deep sequencing of the mRNA and miRNA transcriptomes. Differentially expressed genes were identified using Wald's test, and principal component analysis (PCA) of these genes revealed visible separation of PDACS from NETs and some, but not so clear, second principal component separation of PNETs from INETs ( Figure 3).  (Table S2). There were 1235 common to PNETs and INETs compared with PDACs.
There were 798, 1659, and 54 differentially expressed mRNAs that exhibited a high ability to discriminate between PDACs and PNETs, between PDACs and INETs, and between PNETS and INETs, respectively, with an area under the receiver operating characteristics curve (AUC -ROC) > 0.9 (Table S2); of these, 92, 494, and 16 discriminated between tumor samples in the respective pair comparisons with an AUC -ROC = 1.0.
RNAseq data sets were functionally analyzed by annotation to the Reactome signaling pathway database. Several Reactome pathways were found to be discriminative between PDACs and PNETs. The genes that were significantly downregulated in PNETs compared with PDACs could be classified into 53 functional groups (Table S3). The first 10 significantly overrepresented Reactome pathways included genes from five groups: interleukin-10 signaling, immune system, extracellular matrix organization, degradation of the extracellular matrix, and syndecan interactions (Table S4). Upregulated transcripts could be described with 25 functional groups. Here, the first 10 significantly overrepresented Reactome pathways included genes from five groups: neuronal system; GABA synthesis, release, reuptake and degradation; neurotransmitter receptors and postsynaptic signal transmission; GABA receptor activation; regulation of gene expression in beta cells (Tables S5 and  S6).
On average, 1,161,543 reads mapped to miRBase were obtained per library, representing 41% of total reads. Altogether, 1532 mature miRNAs were detected, of which 452 generated 10 reads on average. The pair-wise comparisons identified in total 258, 297, and 129 significantly dysregulated miRNAs (adjusted p-values < 0.05) between PDACs and PNETs, PDACs and INETS, and PNETs and INETs, respectively (Table S7). There were 157 miRNAs common to PNETs and INETs when  (Table S2). There were 1235 common to PNETs and INETs compared with PDACs.
There were 798, 1659, and 54 differentially expressed mRNAs that exhibited a high ability to discriminate between PDACs and PNETs, between PDACs and INETs, and between PNETS and INETs, respectively, with an area under the receiver operating characteristics curve (AUC − ROC) > 0.9 (Table S2); of these, 92, 494, and 16 discriminated between tumor samples in the respective pair comparisons with an AUC − ROC = 1.0.
RNAseq data sets were functionally analyzed by annotation to the Reactome signaling pathway database. Several Reactome pathways were found to be discriminative between PDACs and PNETs. The genes that were significantly downregulated in PNETs compared with PDACs could be classified into 53 functional groups (Table S3). The first 10 significantly overrepresented Reactome pathways included genes from five groups: interleukin-10 signaling, immune system, extracellular matrix organization, degradation of the extracellular matrix, and syndecan interactions (Table S4). Upregulated transcripts could be described with 25 functional groups. Here, the first 10 significantly overrepresented Reactome pathways included genes from five groups: neuronal system; GABA synthesis, release, reuptake and degradation; neurotransmitter receptors and postsynaptic signal transmission; GABA receptor activation; regulation of gene expression in beta cells (Tables S5 and S6). On average, 1,161,543 reads mapped to miRBase were obtained per library, representing 41% of total reads. Altogether, 1532 mature miRNAs were detected, of which 452 generated 10 reads on average. The pair-wise comparisons identified in total 258, 297, and 129 significantly dysregulated miRNAs (adjusted p-values < 0.05) between PDACs and PNETs, PDACs and INETS, and PNETs and INETs, respectively (Table S7)  There were 126, 250, and 117 differentially expressed miRNAs that were highly discriminative between PDACs and PNETs, between PDACs and INETs, and between PNETS and INETs, respectively, with AUC-ROC 0.9 (Table S8). Of these, 26, 176, and 89 miRNAs discriminated between tumor samples in the respective pair comparisons with an AUC-ROC = 1.0.

Integrated Analysis
The PLS-DA analysis resulted in models with high discriminatory powers (AUC = 1) for both miRNA and mRNA in the second component (Table 4, Figure 5). There were in total of 25 contributors to the first and second components from mRNA and 15 from miRNAs ( Figure 6). mRNAs differentiated mostly PDACs and PNETs from other groups, while miRNAs contributed mostly to PNET and INET differentiation. In addition, 220 correlations with absolute coefficient values above 0.4 were found between contributors to the first three components, of which 175 were negative correlations. This might suggest negative regulation by miRNAs (Table S8, Figure 1).  There were 126, 250, and 117 differentially expressed miRNAs that were highly discriminative between PDACs and PNETs, between PDACs and INETs, and between PNETS and INETs, respectively, with AUC-ROC 0.9 (Table S8). Of these, 26, 176, and 89 miRNAs discriminated between tumor samples in the respective pair comparisons with an AUC-ROC = 1.0.

Integrated Analysis
The PLS-DA analysis resulted in models with high discriminatory powers (AUC = 1) for both miRNA and mRNA in the second component (Table 4, Figure 5). There were in total of 25 contributors to the first and second components from mRNA and 15 from miRNAs ( Figure 6). mRNAs differentiated mostly PDACs and PNETs from other groups, while miRNAs contributed mostly to PNET and INET differentiation. In addition, 220 correlations with absolute coefficient values above 0.4 were found between contributors to the first three components, of which 175 were negative correlations. This might suggest negative regulation by miRNAs (Table S8, Figure 1).

Discussion
PNETs are relatively uncommon neoplasms, constituting 1-5% of pancreatic tumors [9,23,24]. Two major categories of PNETs are well-differentiated NETs and poorly-differentiated neuroendocrine carcinomas (NECs) [25]. In contrast to aggressive PDACs with short survival, PNETs represent neoplasms of a relative indolent clinical course characterized by slow tumor growth [25]. Although a possible common origin between PDACs and PNETs may exist [26][27][28], deep differences between well-differentiated NETs and endocrine carcinomas, as well as adenocarcinomas, have been reported with regard to mutations and gene expression [25,29].

Discussion
PNETs are relatively uncommon neoplasms, constituting 1-5% of pancreatic tumors [9,23,24]. Two major categories of PNETs are well-differentiated NETs and poorly-differentiated neuroendocrine carcinomas (NECs) [25]. In contrast to aggressive PDACs with short survival, PNETs represent neoplasms of a relative indolent clinical course characterized by slow tumor growth [25]. Although a possible common origin between PDACs and PNETs may exist [26][27][28], deep differences between well-differentiated NETs and endocrine carcinomas, as well as adenocarcinomas, have been reported with regard to mutations and gene expression [25,29].
Whole-genome sequencing of sporadic PNETs, mostly of G1-G2 grade, revealed that MEN1, DAXX, and ATRX were mutated in 36.7%, 22.4%, and 10.2% of tumors, respectively, and that PTEN, SETD2, ATM, MTOR, NUMA1, TP53, and KMT2C were mutated at a frequency between 7.1 and 3.1% [21]. Of these genes, MEN1 and DAXX were found to be mutated in one-fifth of our PNET samples; other genes, including ATRX, PGFRA, MTOR, ATM, TSC2, and VHL, were mutated in individual samples. In summary, we demonstrated that the mutational burden of NETs did not differ from that of PDACs, and while some of the mutated driver genes were shared between all three tumor types, others differentiated one tumor type from the other two.
In contrast to genetic analysis, deep sequencing of the mRNA and miRNA transcriptomes identified much more visible separation between PDACS and NETs at the molecular level. Of 20,520 mRNAs and 1532 mature miRNAs identified by RNAseq, 2082, 2621, and 107 mRNAs and 258, 297, and 129 significantly dysregulated miRNAs differentiated PDACs and PNETs, PDACs and INETS, and PNETs and INETs, respectively. In addition, 220 correlations were found, mostly negative, between relative levels of mRNAs and miRNAs for variables that contributed to three components in DIABLO analysis.
The rarity of PNETs has limited the number of large-scale transcriptomic studies, in contrast to numerous studies previously conducted in pancreatic cancer. Bioinformatics analysis of the two data sets from ICGC and two from the NCBI GEO database has identified the five common differentially expressed genes in PNETs (ABCC8, PCSK2, IL13RA2, KLKB1, and PART1) related to their development and prognosis [38]. A microarray-based study, performed on 72 primary PNETs, 7 matched metastases, and 10 normal pancreatic samples, has found that TSC2 and PTEN-two important inhibitors of the Akt-mTOR pathway-are downregulated in PNETs, and FGF13 is upregulated in metastatic compared to nonmetastatic primary tumors [39]. Both mouse and human well-differentiated islet/insulinoma tumors express genes characteristic of mature islet cells, while poorly differentiated tumors associated with liver metastases express genes known to regulate early pancreatic development [40]. Another microarray-based analysis of gene expression patterns in PNETs has shown that upregulated, differentially expressed genes are related to several pathways, including type 2 diabetes mellitus, Ca 2+ signaling, long-term potentiation, and long-term depression, and that downregulated genes are enriched in pancreatic secretion, protein digestion and absorption, and other metabolic pathways. Interferon-stimulated gene protein 15, somatostatin, and synaptosomal-associated protein 25 kDa are identified as hub proteins [41].
When our mRNAseq data sets were annotated according to Reactome signaling, downregulated transcripts in PNETs, when compared with PDACs, were classified into 53 functional groups, including interleukin-10 signaling, immune system, extracellular matrix organization, degradation of the extracellular matrix, and syndecan interactions. Upregulated transcripts fell into 25 groups, including the neuronal system; GABA synthesis, release, reuptake, and degradation; neurotransmitter receptors and postsynaptic signal transmission; GABA receptor activation; regulation of gene expression in beta cells (Tables S4 and S6). These findings proved differences in the biology of well-differentiated neuroendocrine tumors and pancreatic adenocarcinomas, which might relate to different tumorigenesis pathways.
In fact, NETs originate from the diffuse neuroendocrine system, and functional NETs can secrete different peptides and amines categorized as nonspecific markers (e.g., chromogranin A, neuron-specific enolase, and pancreatic polypeptide) and specific markers (e.g., serotonin, gastrin, insulin, glucagon, somatostatin, and vasoactive intestinal peptide). Our functional annotations of differentially expressed genes confirmed significant biological differences between PNETs and PDACs.

Patients
A total of 21 patients were enrolled in this study from January 2018 to December 2018, including 9 with PNETs, 8 with PDACs, and 4 patients with intestinal neuroendocrine tumors (INETs); of the latter, one tissue sample was obtained from liver metastasis. All cases had adequate clinical and pathologic information. The patient's sex, age, tumor site, size, stage, histology, grade, type of surgery performed, additional treatment, and follow-up data were recorded. Tumor grade and stage were evaluated according to the criteria of the 8th edition of the AJCC staging manual (2017) (www.cancerstaging.org ajcc@facs.org). All subjects provided written informed consent prior to participation, and the study complied with the Declaration of Helsinki. PNET patients. The median age of PNET patients (8 women and 1 man) was 45.4 years (range 24-70). The median tumor size was 2.9 cm (range 1.3-3.5); five tumors were located in the pancreatic tail. Two tumors were G1, six were G2, and one was G3, respectively. Six tumors were at stage II, one at stage III, and one at stage IV. In all, but one, patients, distal pancreatic resection or pancreaticoduodenectomy was performed.
PDAC patients. The median age of PDAC patients (6 women and 2 men) was 60.3 years (range . The median tumor size was 2.7 cm (range 1.7-3.3); six tumors were located in the pancreatic head. All tumors were classified as ductal adenocarcinoma. Five tumors were G2, and six tumors were at stage III and 2 at stage IV. In five patients, potentially curative surgery was performed, and in three others, the only laparotomy.
INET patients. The median age of the INET patients (one woman and three men) was 64 years (range 59-67). The median tumor size was 2.6 cm (range 2-3). Three primary tumors were G1, and one liver metastatic tumor was G3. There was one tumor at stage II, one at stage III, and one at stage IV.
All tissue samples for molecular analysis were unfixed postoperative specimens, and all, but one, were primary tumors. The tissue was frozen in liquid nitrogen and stored until use at −70 • C.

Nucleic Acid Isolation
Total nucleic acids were isolated from tumor tissue specimens using the ALLPrep DNA/RNA Micro Kit (Qiagen, Hilden, Germany), following the manufacturer's protocol. DNA sample concentrations were measured using a NanoDrop spectrophotometer, following the manufacturer's instructions, and the DNA was stored at -20 • C. The purity and quantity of total RNA were measured using a Qubit Fluorometer and assessed using an Agilent 2100 Bioanalyzer with an Agilent RNA 6000 Nano Kit (Agilent, Santa Clara, CA, United States). RNA samples were stored at -70 • C.

Profiling of DNA Variants
Ion AmpliSeq Comprehensive Cancer Panel (ThermoFisher Scientific, Waltham, MA, United States) libraries were prepared from DNA samples for the analysis of the coding regions of 409 oncogenes and tumor suppressor genes by sequencing on an Ion Proton sequencer with the Ion PI Hi-Q Chef Kit, loading four samples onto an Ion PI Chip, as described previously [54].

Data Analysis and Variant Calling
Raw reads were processed using the Torrent Suite analysis pipeline version 5.10 and mapped to human genome assembly hg19 using TMAP [88]. Variant calls were made using Ion Reporter Software [89] (version 5.10), using the default parameters for somatic variants. Called variants were filtered using bcftools (version 1.3) and parameters provided in Table S9. Variants were annotated, and their functional consequences were predicted using Cravat [90] version 4.3. Only rare (i.e., with minor allele frequency less than 2%) and non-synonymous variants were analyzed further. From these variants, potential driver mutations were chosen based on one of the following criteria: (i) the CHASM p-value after FDR correction was smaller or equal to 0.1 for a single variant, (ii) the variant was a non-synonymous variant in a known driver gene, or (iii) the variant was pathogenic according to the VEST algorithm.

mRNA Data Analysis
Raw reads were processed using the Torrent Suite analysis pipeline and mapped to the AmpliSeqTranscriptome version of the human genome assembly hg19 by use of TMAP. Reads corresponding to each gene were counted using htseq-count [91] version 0.6. Normalization and differential expression were conducted using DESeq2 [92] version 1.18, using the default parameters and options. Differentially expressed genes in the pair-wise comparisons were identified using Wald's test. The resulting p-values were adjusted for the testing of multiple hypotheses using the Benjamini-Hochberg procedure to control the false discovery rate. The false discovery rate threshold was set to 0.05, and only probe sets exhibiting a minimum 2-fold change in mean relative expression were included in the functional analysis. Overrepresentation of gene ontology and Reactome pathways with hypergeometric testing was determined using ClueGO version 2.5.4, with false discovery rate p-value correction [93].

Integrated Data Analysis
Transcriptome and miRNA data integration (including PLS-DA and correlation computation) were conducted using DIABLO with mixOmics functions [97] package (mixOmics version 6.8.5). The tuning of the component and the variable number was performed according to the tutorial at http://mixomics.org/mixdiablo/case-study-tcga/.

Conclusions
To sum up, our study was the first to compare directly and simultaneously mutational and transcriptional profiles of well-differentiated neuroendocrine tumors and ductal adenocarcinomas of the pancreas using massively parallel sequencing. We confirmed the presence of two driver genes, MEN1 and DAXX, mutations in which were found in one-fifth of PNET samples but not in other tumor types. Among the most-mutated genes, KRAS and TP53 were mutated only in cancers, while five other genes (HNF1A, CREBBP, RNF213, CASC5, and APC) were mutated with similar frequencies in adenocarcinomas and neuroendocrine tumors. However, most other driving mutations were found only in individual tumor samples or in small percentages of tumors. In contrast to the genetic findings, the clearest differences were transcriptomic. While differentially expressed genes may indicate different sites of origin for adenocarcinomas and neuroendocrine tumors, they cannot be linked to significant differences in the aggressiveness of PNETs and PDACs. Finally, our study revealed that both mRNA and miRNA profiles had a high ability to discriminate between tumor samples. In addition, integrative analysis using the novel algorithm DIABLO allowed us to construct a multi-omic biomarker model, based on quantitative data from miRNA and mRNA. Both profiles showed high discriminatory powers, mainly in the second component, where they reached AUC of 1 for all the groups. Unsurprisingly, the miRNA profile was also highly correlated to mRNA expression. Hence, miRNAs can be considered as a new class of pancreatic tumor biomarkers due to their high stability.

Conflicts of Interest:
The authors declare no conflict of interest.