Genetic datasets of PNETs using whole exome sequencing (WES), whole genome sequencing (WGS), and targeted gene-panels of more than 400 genes were searched for in Pubmed on 25^th of September 2023 using the following search criteria: “Pancreatic neuroendocrine tumor” AND (“sequencing” OR “NGS”), “Pancreatic neuroendocrine neoplasm” AND (“sequencing” OR “NGS”),”Insulinoma” AND (“sequencing” OR “NGS”), “Glucagonoma” AND (“sequencing” OR “NGS”),”Somatostatinoma” AND (“sequencing” OR “NGS),”Gastrinoma” AND (“sequencing” OR “NGS”),”VIPoma” AND (“sequencing” OR “NGS”), “Serotonin-producing tumors” AND (“sequencing” OR “NGS”), “ACTH-producing tumors” AND (“sequencing” OR “NGS”), “Gastroenteropancreatic” AND (“sequencing” OR “NGS”), “GEP-NET” AND (“sequencing” OR “NGS”), “GEPNET” AND (“sequencing” OR “NGS”), “PNET” AND (“sequencing” OR “NGS”), “P-NET” AND (“sequencing” OR “NGS”), “PanNET” AND (“sequencing” OR “NGS”), “Pan-NET” AND (“sequencing” OR “NGS”), “PNEN” AND (“sequencing” OR “NGS”). The studies were screened by one reviewer and included datasets were also retrieved by the same reviewer.

All pancreatic functional tumor types in the search were retrieved from Guilmette et al (2019) (24) and WHO Classification of Tumors of Endocrine Organs (2017) (25). The search included words appearing in the title and abstract of each article. Selected articles were then screened individually for inclusion and exclusion criteria as stated below.

We included studies using WES, WGS and/or a gene panel of more than 400 genes to determine the somatic variants in PNETs. The cut-off of 400 genes was chosen based on a Japanese study with use of data from The Cancer Genome Atlas (26), showing that cancer-panel sizes at or above 400 genes had the best statistic performance power to identify tissue mutations in comparison to WES (26). In our dataset, we included human PNETs with tumor grade G1 or G2. Wd-G3 PNETs were not included in the present analyses given the interpretation difficulties and the fact that G3-NETs were not separated from G3 neuroendocrine carcinomas (NECs) in articles based on WHO classifications prior to 2017 (13). Both patients with and without tumor syndromes with both functioning and non-functioning PNETs were included. Both studies using formalin-fixed paraffin-embedded (FFPE) and fresh frozen samples were included.

Exclusion criteria included studies on cell-lines and human studies on extra-pancreatic tissues, carcinomas, non-islet PNETs, metastatic tissues, PNETs with tumor grade G3 only and NECs; Studies with data without proper sample characteristics were excluded. Studies without obtainable full text or main dataset were also excluded. Replicated data were only presented once from the original study. The literature search flow chart is presented in Figure 1 . The full list of excluded and included papers is listed in Supplementary Tables 1 , 2 and Table 1 , respectively.

The meta-analysis of included data involved the following steps: 1) Select all datasets in each published database comprising samples belonging to the chosen conditions. 2) For each separate dataset, find and select PNETs with a tumor grade of G1 or G2. 3) Categorize the samples based on tumor type (functional vs non-functional) and sporadic vs hereditary when specified, retrieve DNA variations, including mutations, and copy number variations (CNVs; including amplifications (Amp), deletions (Del), and copy-neutral loss of heterozygosity (cnLOH) of gene region or whole chromosome), where applicable, 5) exclude synonymous variants, 6) exclude intronic variants, 7) exclude common variants identified in the general population (GnomAD, frequency < 1%), 8) list the mutations and CNVs in genes somatically altered in PNETs, 9) compose data to a single mutation annotation format MAF-like file for analysis as described below.

One dataset was aligned to the reference genome GRCh38 and was converted to GRCh37 using LiftOver (https://genome.ucsc.edu/cgi-bin/hgLiftOver). Twelve intronic variants in CUX could not be converted. Alterations on nucleotide level were predicted using the software Transvar (37) using the amino acid alteration and coding sequence (CDS) as input. The generated nucleotide alteration (reference and alternative allele) was compared to four of the included datasets (14, 15, 19, 32), n=71/225. By using the software, both strands could be considered altered, thus we accepted the generated nucleotide if it followed the base pairing rule in comparison to the datasets. Two studies lacked HGVS description of amino acid alterations (32, 33), and these were generated using the genomic coordinates of the CDS and the represented amino acids (reference and alternative amino acid) or nucleotide change, respectively. The study of Scarpa et al. (13) used a different annotation when describing indels, and the HGVS annotation was achieved using Transvar based on the genomic coordinates of the CDS and the chromosome number. A few splice-variants were excluded due to missing data.

Data were collected from each included article and the data were listed in a MAF-like structure. The collection of data was achieved in excel and the data-frame was compiled to a CSV-file before the generation of the MAF-file in R. The resultant MAF-like file was generated for all samples and consisted of the following information: Genetic data, including Human Genome Organization (HUGO) gene nomenclature symbol, variant classification, variant type, reference allele, tumor seq allele1, tumor seq allele2, and amino acid change along with the tumor sample barcode and the respective article. Clinical data about metastasis, tumor type, and tumor syndrome was also added to the MAF file for each tumor sample barcode. When applicable, CNV data were collected and added to an additional customized table in a second file. The file could be merged to the MAF-like file during the meta-analysis.

The generated MAF-like file including all retrieved sample information was used for further analysis. We categorized each PNET of a patient into sporadic vs. hereditary, including germline variations in MEN1, VHL, NF1, CDKN1B, BRCA2, MUTYH, CHEK2, PTEN, by use of the descriptions from the articles and corresponding datasets. Further, each PNET was categorized into non-functional (e.g. PPomas, and other PNETs without hypersecreting syndrome) vs. functional PNETs (e.g. insulinoma, glucagonoma, VIPoma, gastrinoma) according to the retrieved published data. Further, the presence of distant metastases was registered, if published.

The R-package Maftools (38) was used to produce a summary of the file, oncoplots and a plot of implicated oncogenic pathways, along with analyses of gene-specific mutations, co-occurrence of mutations, and variant enrichment.

The pathogenicity of missense mutations in DAXX and ATRX were assessed using Polyphen-2 (39). Protein domains in DAXX and ATRX were retrieved from Wang et al. (40).

P-values were determined using two-tailed Fishers-exact test, a Freeman-Halton extended Fishers exact test, or Mann-Whitney U test where appropriate. P-values less than 0.05 were considered significant. Data were analyzed using R (v. 4.2.3.).

Our study included 14 of 300 datasets available using our search criteria, Table 1 . The 14 datasets comprised a total of 225 individual G1/G2 PNET samples from 221 patients. Of the 225 PNETs included, 162 (72%) were sporadic with no associated germline alteration and 30 (13%) were hereditary with germline predisposing alterations in MEN1 (14 PNETs in 11 patients), MUTYH, CHEK2, VHL, CDKN1B, BRCA2, and PTEN, Table 2 and Supplementary Table 3 . Two patients harbored variants in two predisposing genes (MEN1/MUTYH and CHEK2/MUTYH). In 33 PNETs (15%), information on germline predisposition was missing.

Patients with germline predisposition were significantly younger than patients with sporadic PNETs, 46 vs. 56 years, p=0.005. The PNETs were functional in 45 (20%), of which insulinomas accounted for 37 (82%), non-functional in 156 (69%), and with unspecified functionality in 24 (11%). Our meta-analysis showed no difference in functional and non-functional PNETs in hereditary vs. sporadic PNETs, 5/22 (22.7%) vs. 40/146 (27.4%), p=0.80.

Of the 225 PNETs, 25 (11%) were accompanied with distant metastasis. No difference was seen between genetically predisposed and sporadic PNETs regarding the presence of metastasis.

The 225 PNETs had a total of 6,194 reported non-synonymous somatic variants. We identified a difference between computed nucleotide change by Transvar and the annotated change in 74 cases, accounting for 7.3% of tested nucleotide alterations. In addition, the Transvar software failed to identify nucleotide change for 420 of the reported changes. In 118 cases, the cause of failing was an invalid gene annotation, 294 lacked a valid transcript, 6 had an invalid mutation nomenclature and 2 were out of range. Due to the pronounced difference between the Transvar output and the published data, we decided only to use the annotations described by the articles regarding the nucleotides on each tumor allele. The overall distribution of transversions and transitions in the resultant dataset is shown in Supplementary Figure 1 .

The most frequently somatically altered genes in all PNETs are highlighted in Figure 2 . The top three variated genes were MEN1 (n=95; 42%), DAXX (n= 37, 17%) and ATRX (n=27; 12%). Somatic point mutations/indels in MEN1 were identified in 67 PNETs and CNVs in 41 PNETs, including 13 PNETs with both a point mutation/indel and CNV in the gene. The somatic MEN1 point mutations/indels in the coding region were distributed throughout the gene as shown in Figure 3 (lower row).

Variants in the coding regions of DAXX and ATRX are presented in Figure 4 . The DAXX and ATRX mutations were unique except for one DAXX alteration, which was present in two PNETs. All of the missense mutations in DAXX (n=3) and ATRX (n=6) were predicted as probably damaging by PolyPhen2. Notably, the DAXX missense mutations were located in close proximity to each other (amino acid 328 to 331) in the histone-binding domain. The majority of missense variants in ATRX (5/6) were all located in the ATPase domain spanning residues 1,550–2,226.

MEN1 point mutations was identified in 30% of the analyzed tumors and 42% when including CNVs. Consistently, previous smaller individual studies with or without CNV analyses reported a frequency of 25–44% (11–14). In PNETs from MEN1 patients, MEN1 was somatically altered as a second hit in all 13 with available data. Inactivation of a second MEN1 allele may happen due to LOH, point mutations, or epigenetic inactivation. In our dataset, 92% (12/13) of the MEN1 patients acquired LOH in the MEN1 locus as the somatic second hit (one patient had no MEN1 somatic data). Likewise, LOH has earlier been proposed to be the main mechanism of full MEN1 inactivation, accounting for 90% of PNETs in MEN1 patients (41). MEN1 patients did not share other somatic gene mutations, which may reflect the sovereignty of menin deficiency in their PNET development.

In sporadic PNETs, Jiao et al. identified point mutations/indels in MEN1 in 30/68 (44%) of the examined cohort. (12). In our meta-analyses, which did not include this study, 52/162 (32%) sporadic PNETs harbored somatic MEN1 point mutations/indels. Only six sporadic PNETs had a CNV of the MEN1 locus, and CNVs in MEN1 were significantly more frequent in hereditary PNETs compared to sporadic PNETs. This difference may, however, be assigned the low number of CNV data in sporadic PNETs. Somatic MEN1 alterations were equally distributed in non-functioning vs. functioning PNETs with a high frequency of CNVs in each group. This further highlights the need of CNV analyses in all subtypes of PNETs when studying their tumorigenesis.

DAXX and ATRX were the second and third most frequently altered genes in our meta-analysis. Two patients with hereditary MEN1 had a somatic alteration in both MEN1 and DAXX. The heterogeneity of ATRX and DAXX mutations has been explored earlier and it was also apparent in the cohort of PNETs. Only one variant in DAXX was altered in two PNETs, the remaining DAXX-mutations were exclusive and specific to each PNET.

A mutational hotspot in DAXX has earlier been proposed, resulting in alterations in the amino acid NM_001141970.1, Serine 102; p.(Ser102Leu) and p.(Ser102*) (42). In our dataset, which did not include the mentioned study, we did not identify variants altering Serine 102. Instead, we identified another hotspot for three non-synonymous missense variations in DAXX spanning from amino acid 328 to 331, located in a histone binding domain (amino acids 178–389) of the protein (40). Indeed, Jiao et al. (12)also identified a DAXX missense variant in this domain p.(Ala297Pro) and an in-frame deletion p.(Leu309GlnPhe), and Quevedo et al.(43 ) identified a p.(Ala313Val) variant in a PNET. Taken together, the DAXX histone-binding domain 178–389 represents a new mutational hotspot in PNETs, although DAXX mutations were widespread.

ATRX mutations were also unique in regard to the position in the gene and were only observed once and with no tendency toward specific mutation types. The six somatic missense mutations in ATRX were all predicted probably damaging. Five of the identified missense mutations, (p.(His1759Asp), p.(Met1800Ile), p.(Lys2036Glu), p.(Ser2116Phe), and p.(Tyr2176Cys), were clustered in the C-terminal ATPase domain (amino acid 1,550–2,226) of ATRX (40). This domain contains seven conserved motifs responsible for the ATP-hydrolysis and may further be responsible for DNA translocase activity (44, 45). Experimental evidence has shown that this domain is DNA-dependent and may be associated with chromatin remodeling or DNA replication (46). The variant p.(Lys2036Glu) is placed within the highly conserved helicase domain IV (47, 48). The five mutations could have a noticeable effect on the function of ATRX, and may hinder the ATPase activity of the protein. The mutational hotspot may indicate, that this function is pivotal for endocrine cells, as the interaction with DAXX at H3.3 still may be intact. Interestingly, other PNET study data from Jiao et al. and Quevedo et al. showed a total of four ATRX missense variants, all clustering in the ATPase domain (12, 43), supporting the importance of the ATRX ATPase activity in PNETs.

In the sporadic cohort, DAXX and ATRX were mutually exclusive except for one PNET which had a CNV in both genes. This further strengthens the potential for DAXX and ATRX alterations as independent tumor drivers in PNETs. Furthermore, DAXX was altered significantly more in non-functioning compared to functioning PNETs, indicating a frequent role of DAXX in the tumorigenesis of non-functioning PNETs. In functioning PNETs, only one DAXX and one ATRX alteration was seen in a glucagonoma, and an insulinoma, respectively. ATRX and DAXX mutations has previously been associated with a more alpha cell-like phenotype in a study not included in the present meta-analysis (49).

YY1 mutations was especially frequent in insulinomas, in which 11% harbored a point mutation or indel. Mutations in YY1 was less frequently seen in the meta-analysis compared to the included Asian study of Cao et al., where up to 30% of insulinomas had YY1 alterations (15). The recurrent YY1 mutation p.(Thr372Arg) occurred in 14% in our analysis, narrowing the frequency of 8–30% of reported in three individual studies, of which only the study by Cao et al. was included in our meta-analysis (15–17). Of note, Lichtenhauer et al. found a lower prevalence of YY1 alterations of 12% in insulinomas from their Caucasian cohort (16), compared to sporadic insulinomas from an Asian cohort (15).

Surprisingly, YY1 was not exclusively altered in insulinomas in our datasets. One non-functioning tumor also harbored the characteristic p.(Thr372Arg) mutation, and the YY1 region was amplified in nine, especially non-functioning, PNETs. The latter may indicate that YY1 could be a proto-oncogene just as observed for the gene amplifications and overexpression of MDM2 in cancer (50). YY1 codes for the potential proto-oncogenic ying yang 1 transcription factor (51), confirming the relevance of YY1 in tumorigenesis.

In hereditary PNETs from other than MEN1 patients, germline mutations were observed in MUTYH, CHEK2, BRCA2, VHL, PTEN, and CDKN1B, of which one had germline mutations in both CHEK2 and MUTYH. Furthermore, one MEN1 patient had germline mutations in both MEN1 and MUTYH. Moreover, DAXX variants occurred in four patients with CHEK2 and/or MUTYH germline presentation, highlighting the need of expanded genetics in PNET tumorigenesis research.

Of the six patients with germline MUTYH variants, three of the variants were predicted as pathogenic (13). Our datasets indicated a high level of tumor mutational burden in these individuals, even though the variant numbers between datasets should be compared with caution.

While the MUTYH protein is involved in base-excision-repair during DNA damage (52), BRCA2 and CHK2 are members of the homologous recombination pathway, repairing double-stranded DNA breaks (53). A BRCA germline variant was detected in one patient, and CHEK2 germline variants were identified in four patients. Genomic instability caused by defective DNA repair proteins is a well-described hallmark of cancer (54). Accordingly, patients with germline variants in MUTYH, BRCA2 or CHEK2 should be considered at risk of PNET development, as well as at risk for hereditary cancer (55). The presence of BRCA2 germline variants in pancreatic lesions, including PNETs, has been described earlier (13, 14, 56–58).

Germline mutations in VHL causing von-Hippel Lindau syndrome appeared in five patients from our datasets. Development of PNETs in this rare syndrome is seen in 17% (273/2,330) of all VHL patients according to The European-American-Asian-VHL-panNET-Registry (59, 60). Especially variants in exon 3 have been associated with malignancy (59, 61). Of the five patients with a VHL germline alteration, one had a missense variant in exon 3 and somatic LOH of the second allele and was metastasis-free (13). Somatically, we identified eight VHL variants in the whole dataset, of which one was a frameshift mutation and seven were missense variants distributed in all three exons. However, only a single of these PNETs had metastasis and this variant was in exon 2. Seven of the eight PNETs with VHL variants were non-functioning. This predominance of non-functioning PNETs is in line with other studies (61, 62).

Notably, one PNET-patient included in our datasets, had been diagnosed with Cowden syndrome, harboring a deleterious germline PTEN-variant (34). A case report of the patient has been published (63), highlighting the rarity of the association of PNET development and the syndrome. Lastly, one PNET harbored a germline variant in CDKN1B. Loss of germline CDKN1B is referred to as MEN4 syndrome and can cause a MEN1-like phenotype. Only very few PNETs have been described having germline variants in CDKN1B (64, 65). Notably, patients with hereditary PNETs were significantly younger than patients with sporadic PNETs, as described earlier (1, 66).

PNETs, which are mostly non-functioning, have been associated with variants in AKT/PI3K/mTOR pathway, i.e. PTEN, TSC1, TSC2, PIK3CA and DEPDC5, which regulates cell survival and proliferation (13, 67). It has been suggested that patients with mTOR pathway aberrations may benefit from therapy using mTOR inhibitors (12, 13). In our combined datasets, we also found variants in the PI3K/mTOR pathway, especially in the non-functioning tumors (30/156 PNETs). The most abundant altered genes were coding for PTEN (n=9), TSC2 (n=6), TSC1 (n=5) and MTOR (n=4). Only one tumor had mutations in more than one gene encoding a PI3K/mTOR pathway interactor. Interestingly, the use of everolimus in non-functioning NETs improved the progression-free survival compared to placebo (11 vs 3.9 months, RADIANT-4 trial, (68)), and although the documented use of everolimus in functioning NETs is sparse, mTOR-inhibitor treatment may also affect the clinical symptoms, independent of restrain of tumor growth (69, 70).

RTK-Ras (n=20), Wnt pathway (n=16) and NOTCH (n=16) pathways were also altered in more than 10% of non-functioning PNET. While most of the NOTCH associated proteins were tumor suppressors, the RTK-Ras pathway mostly consists of proto-oncogenes ( Supplementary Figure 2 ), and a maximum of three non-functioning PNETs had mutations in the same gene (ERBB2, AXIN2, and NOTCH2), thereby indicating common pathway variations, but not necessarily gene specific alterations also among the non-functioning tumors. The identification of aberrant pathways may in turn give insight into potential targets for drug designs in the future. Drugs targeting the distinct protooncogenic pathways RTK-Ras, Wnt and NOTCH do already exist. These include for example the multi-kinase inhibitor sorafenib (71), which may impact the entire pathway. For the Wnt and NOTCH-pathways only very few FDA approved drugs are available, however a few new pathway-specific drugs have entered clinical trials (72, 73).

Strengths our study included the systematic search for eligible articles and, by nature of a meta-analysis, our large sample size from 14 datasets compared to individual datasets. Moreover, our inclusion criteria of reported G1/G2 PNETs ensured a focus on wd-PNETs, as G3 PNET is a newly defined entity, not incorporated in most previous studies.

In the WHO version of 2010, wd-NETs of low to intermediate grade was defined as G1-G2, whereas poorly differentiated, high-grade neoplasms were defined as G3 (74). In 2017, WHO altered the definition and introduced a G3 NET described as well-defined with a distinct morphology compared to G3 neuroendocrine carcinomas (NECs) (75). G3-NET shares a molecular mutation profile similar to G1 and G2 NET and are associated with better survival than the poorly differentiated G3 NEC (76), which also has a mutational profile more similar to pancreatic ductal adenocarcinomas (77). G3 NETs may be difficult to distinguish from (G3) NECs, as both entities share a high proliferation rate (Ki67 proliferation index > 20% and/or mitotic rate above 20 per high power field) (78). New studies on wd-PNETs should, however, include G3 PNETs according to the new WHO classification.

Limitations included the varying methods in the included articles, which may have affected the results. The variations include different tissue preparations (FFPE or frozen), different type of library preparations (with or without gene enrichment), different sequencing platform, and resultant different read depth. Six of the articles did not describe the average read depth, although for one article the quality of the sequencing was described using Q₃₀ (27), which indicates the likelihood of incorrect base calling in 1 of 1000 times (79). The differences between the article’s methods may be reflected by the variant calling, as some articles presented with very few variants (28, 29, 34) compared to others (32, 33, 80). Of note, after retrieving the variants, we excluded synonymous and intronic variants, even though such variants can alter the splicing of genes (81).

Moreover, CNV information could only be gathered from 33% of the PNETs, of which we only included structural variations, i.e. amp, del, and cnLOH as structural variations, e.g. tandem duplications, inversions, and translocations, have unpredictable consequences and such variations were only described in one of the 14 included studies (n=33) (32). Underreporting of clinically important CNVs could, therefore, not be ruled out.

Lastly, patients without a known germline predisposition were termed ‘sporadic’, even though most articles did not specify other germline gene results but for MEN1. This may lead to underreporting of rarer germline gene variants in PNETs such as VHL, PTEN, CDKN1B, BRCA2, CHEK2, and MUTYH.