In-depth analysis of the genomic landscape of 86 metastatic neuroendocrine neoplasms reveals subtype-heterogeneity and potential therapeutic targets


 Metastatic neuroendocrine neoplasms (mNEN) form clinically and genetically heterogeneous malignancies, characterized by distinct prognoses based upon primary tumor localization, functionality, grade, proliferation index and diverse outcomes to treatment. Here, we report the mutational landscape of 86 whole-genome sequenced mNEN. This landscape revealed distinct genomic subpopulations of mNEN based on primary localization and differentiation grade; we observed relatively high tumor mutational burdens (TMB) in neuroendocrine carcinoma (5.45 somatic mutations per megabase) with TP53, KRAS, RB1, MYC and APC as major drivers versus an overall low TMB in neuroendocrine tumors (1.08). Furthermore, we observed distinct drivers which were enriched with somatic aberrations in pancreatic (MEN1, ATRX, DAXX, PCNT and SETD2) and midgut-derived neuroendocrine tumors (CDKN1B). Finally, 49% of mNEN patients revealed extensions of their treatment-repertoire based upon actionable (and responsive) somatic aberrations; potentially directing improvements in mNEN treatment strategies.


Introduction 50
Neuroendocrine neoplasms (NEN) is a heterogeneous and uncommon tumor type. It can 51 arise from any of the neuroendocrine cells distributed widely throughout the body. 52 Clinically, a distinction is made between the poorly differentiated neuroendocrine 53 carcinomas (NEC) and the better differentiated neuroendocrine tumors (NET) 1 , the latter 54 are further subdivided based on their primary site in pancreas (pNET), gastro-intestinal tract 55 or lung NET. Further distinctions are made based upon grade (as assessed by Ki-67 or MIB-1 56 staining as a measure of proliferation index), differentiation, histology (small-cell vs. large-57 cell) and functionality (the presence or absence of hormone secretion resulting in typical 58 clinical syndromes dependent upon the predominant hormone that is secreted). Tumor 59 grade and differentiation are associated with prognosis, and all the aforementioned factors 60 affect the choice of treatment. However, also in small subgroups of NEN, such as well-61 differentiated low-proliferating pNET, marked clinical and genetic heterogeneity occur, as 62 well as, vastly different responses to treatment. Thus, the parameters by which NEN are 63 currently classified do not sufficiently separate patients and tumors according to prognosis 64 and response to therapy. Nonetheless, certain anti-tumor therapies (i.e., sunitinib and 65 everolimus) have been registered for distinct NEN-subtypes. Hence, there is a high unmet 66 need to better classify and understand these diverse tumors, ultimately leading to more 67 tumor-or patient-tailored therapeutic strategies. 68 69 Thus far, limited whole genome sequencing data are available for NEN, probably reflecting 70 the rarity of this disease. Currently, pNET have been characterized most extensively; 81 71 primary tumors were subjected to whole genome sequencing as part of the PCAWG project 2 72 and another 104 primary pNET were described by Scarpa et al. 3 . Additionally, smaller series 73 have been published containing lung and gastrointestinal NET. 4,5 These studies have shown 74 that NET have a relatively stable genome and only few driver mutations. However, these 75 studies were all performed on primary tumor specimens, whilst a patient generally dies 76 from the consequences of metastatic disease. Additionally, we know from other tumor 77 types that marked heterogeneity can occur between primary and metastatic tumor cells 6-9 , 78 due to inherent genomic instability and/or the influence of targeted or cytotoxic treatment 79 on the tumor genome. These discrepancies should be taken into account when assessing a 80 patient's prognosis and possible treatment options, and can be better understood through 81 P a g e 3 | 35 thorough genomic characterization of metastases. To date, analysis of metastatic NET is 82 limited to two studies describing series of five patients with NET originating in the pancreas 83 and the small intestine (or midgut), respectively. 10,11 These studies have shown focal 84 amplification of MYCN concomitant with loss of APC and TP53 in one sample as important 85 metastatic genetic aberrations. For NECs, only two series of whole genome sequencing of 86 the primary tumors of 1) five cervical and 2) 12 genitourinary NECs have been published. 12,13 87 88 Whole genome sequencing was performed on 86 NEN metastases (mNEN), one per distinct 89 patient, which were biopsied as part of the CPCT-02 study. 14 We report here on the 90 presence of genomic alterations, mutational and rearrangement signatures for the whole 91 mNEN cohort and reveal genomic characteristics and alterations distinguishing mNEC from 92 mNET. Furthermore, we make a genomic distinction between pancreas-and midgut-derived 93 mNET. Additionally, we investigated the presence of actionable genetic alteration within 94 mNEN patients, which might render them eligible for off-label or experimental systemic 95 treatments to extend therapy options. 96

Results 97
Overview of included patients within the CPCT-02 mNEN cohort.

98
A total of 109 patients, originally classified as having a neuroendocrine neoplasm, were 99 included in the CPCT-02 study and had a metastatic biopsy taken in parallel with a blood 100 control (Figure 1). Five patients were excluded because of missing or withdrawn informed 101 consent, and another five had non-evaluable biopsies due to low (<20%) tumor cell 102 percentage or low DNA yield. Thirteen biopsies were excluded because of incomplete 103 clinical records, misclassifications of the tumor (based on additional checks of the medical 104 records), or were duplicate biopsies from the same patient. An overview of the mNEN 105 patient inclusion per participating Dutch center (n = 13) can be found in Supplementary 106 The metastatic tumor biopsies and corresponding peripheral blood controls from the 109 remaining 86 distinct patients were whole genome sequenced using paired-end protocols, 110 to a median mean read coverage of 107x (Q(uartile) 1 -Q 3 : 99x-115x) and 38x (Q 1 -Q 3 : 35-42x), 111 P a g e 4 | 35 respectively to a median in silico estimated tumor cell purity of 0.7 (Q 1 -Q 3 : 0.5-0.81) 112 (Supplementary figure 2a-b). 113

121
The mutational landscape of metastatic neuro-endocrine neoplasms reveals differences related to primary 122 localization and degree of differentiation.

161
High TMB (≥10) are often associated with DNA repair deficiency and/or tumors with 162 sensitivity for immune therapy, e.g. checkpoint inhibitors. Four mNEC samples, all from 163 unknown origin, and a single pancreatic mNET showed this high-TMB genotype (Figure 2a). 164 One mNET displayed signs of BRCA2-associated homologous recombination deficiency 165 (HRD), as determined using the CHORD classifier which is mainly based on deletions with 166 flanking microhomology and 1-100kb structural duplications (Figure 2j; supplementary 167 figure 4). Further inspection revealed that this mNET harbored a somatic frameshift 168 mutation within RAD51C, a known HRD-associated gene. 16 169 Regional hypermutation (kataegis).

170
Regional hypermutation (kataegis) was detected in six mNEC; and SBS4, which is associated with smoking; likely due to tobacco mutagens. This could 211 reflect that these metastases could be primary lung NEC, as a relation between mNEC on 212 non-pulmonary origin and smoking is not known. However, as no somatic coding mutations 213 in canonical lung cancer-associated genes were observed and the clinicopathological data of 214 these patients did not point to any different primary tumor other than a NEC, it seems 215 unlikely that these could be primary non-small cell lung cancers. 216 217 Strikingly, the only high-TMB (pancreatic) NET was strongly characterized by SBS11 which 218 exhibits a mutational pattern resembling that of alkylating agents, with a strong enrichment 219 for C/T (G>A) transitions. Previously, an association between treatment with the alkylating 220 agent temozolomide and SBS11 mutations has been found. 21,23 This same patient showed 221 the highest tumor mutational burden with a TMB of 21.4 (median TMB of NET: 1.1) and was 222 treated with a combination of 5-fluorouracil and streptozocin before undergoing a biopsy 223 for the CPCT-02 study. Streptozocin is a capable of DNA alkylation and inhibition of DNA 224 synthesis, and its mechanism of action closely resembles that of temozolomide. 225

226
One mNET was strongly characterized by SBS36, associated with base excision repair (BER) 227 deficiency due to MUTYH alterations, C>A mutations and previously also seen in pancreatic 228 NET. 23-25 Strikingly, this tumor did not harbor specific somatic alterations within MUTYH but 229 possessed a heterozygous germline pathogenic missense mutation within MUTYH 230 (c.527A>G / p.Tyr176Cys; rs34612342) coupled with a complete loss of a single chromosome 231 1, resulting in subsequent loss-of-heterozygosity. We observed an overall heterogeneous pattern of putative drivers, the most frequently 259 putative driver was found to be CDKN2A/B (n = 18; 15), followed by TP53 (n = 17), CDKN1B 260 (n = 12), PTPRD (n = 11), KRAS (n = 11), MEN1 (n = 11) and RB1 (n = 11). Strikingly, a 261 significant portion of the total mNEN cohort had no mutual putative driver(s) (9 out of 86; 262 10%) and only contained patient-specific putative drivers. 263 We next investigated whether any form of mutational enrichment, such as somatic 265 alterations within certain genes (mutations and/or copy-number alterations) or evidence of 266 large-scale events (kataegis and chromothripsis), could be related to one of our three major 267 subgroups relating to subtype or primary localization; being mNEC (n = 16), pancreas-(n = 268 20) and midgut-derived mNET (n = 40). Using a one-sided Fisher's Exact Test (with 269 Benjamini-Hochberg correction), we detected the enrichment of at least one such event(s) Midgut; Q 1 -Q 3 : 0.75 -1.39) and 1.07 (mNET -Unknown; Q 1 -Q 3 : 0.84 -1.53) to 1.27 (mNET -283 Other; Q 1 -Q 3 : 1.10 -1.44) and 1.35 (mNET -Pancreas; Q 1 -Q 3 : 0.9 -2.12). A similar pattern 284 was detected regarding the number of distinct genes with coding mutations. Midgut-285 derived mNET also presented a surprisingly low number of structural variants compared to 286 the other mNET sub-populations. 287 288 Next, we investigated possible differences in putative drivers between our major mNET sub- investigating the statistically significant large-scale copy-number alterations of the 299 chromosomal arms, we also detect striking differences between the major subgroups 300 (supplementary figure 10). Within mNEC, we detected a large number of samples (69%) 301 harboring a loss of 22 q . Midgut-derived mNET revealed amplifications of chromosome 4 p/q , 302 5 p/q , 7 p/q , 10 p/q , 14 p/q , 20 p/q and loss of 9 p/q in various samples (~30%) and a loss of 18 p/q in 303 66% of samples. This re-confirms the high frequency of chromosome 18 loss in midgut-304 derived NET and the association with DDC 26 , as DCC is the most recurrently mutated gene 305 on chromosome 18 in our cohort also (n = 5). Finally, over half of pancreas-derived mNET 306 revealed amplifications of chromosome 5 p/q , 7 p/q , 9 q , 12 p/q , 13 q , 14 p/q , 17 p/q , 18 p/q , 19 p/q , 20 p/q 307 and loss of 22 q . 308 309 Unbiased driver gene analysis (dN/dS) on midgut-derived mNET presented CDKN1B whilst DAXX, SETD2, CREBBP, PCNT, KDR and TSC2 were found to be mutated only within pancreas-314 derived mNET. Moreover, MEN1 was found to be mutationally enriched within pancreas-315 derived mNET when compared to the entire mNEN cohort (q ≤ 0.05). Several midgut-316 derived mNET (n = 9; 23%) did not readily present a shared mutual driver and only harbored 317 somatic mutations in private or as-of-yet unassociated cancer driver genes. 318 Clinically-actionable mutations.

319
We observed forty-two mNEN (49%) harboring one or more target-specific or general 320 somatic aberrations which are known as possible (and responsive) druggable targets against 321 currently-available (or under development) treatment agents are available. Twenty-one 322 mNEN (24%) harbored somatic aberrations corresponding to a treatment that is currently 323 registered for NEN or specifically for the NEN-subtype of that particular patient (figure 5, 324 supplementary table 1). In addition, twelve patients (14%) could benefit from therapies that 325 P a g e 11 | 35 are off-label, but are commonly considered best practice for NEN. Another eight patients 326 (9%) could benefit from drugs which are registered for another indication but not currently 327 administered in NEN treatment. Additionally, six tumors harbored an aberration rendering 328 them sensitive to a drug that is still in development; including a single patient with no 329 actionable alterations otherwise. We found RB1 (n = 11), KRAS (n = 11), MTAP (n = 5), high-330 TMB (≥10; n = 5), RICTOR (n = 4) and TP53 (n = 4) to be the most frequently observed 331 (target-specific or general) somatic aberrations which granted eligibility to various possible 332 treatment options. In total, ten midgut-derived mNET (25%) and eleven pancreas-derived 333 mNET (55%) revealed potentially responsive alterations in various genes and most strikingly, 334 almost all mNEN (94%) revealed potential responsive targets due to RB1 and/or KRAS 335 mutations or towards checkpoint inhibitors due to high tumor mutational burden (≥10). 336

Discussion 337
Historically, NEN has long been known as a difficult malignancy to diagnose, monitor and 338 treat due to presentation of an inherently wide spectrum of disease progression, cellular 339 differentiation and low mutational burden, resulting in few targetable mutations and a 340 relatively stable tumor genome. Indeed, mNET is characterized by the lowest TMB of all 341 metastatic cohorts sequenced within the CPCT-02 study. 14 This study is the first to have an The single high-TMB pancreas-derived mNET presented a striking contribution of the 355 mutational signature associated with alkylating agents (temozolomide) and was previously 356 P a g e 12 | 35 treated with a combination of 5-fluorouracil and the alkylating antineoplastic agent 357 streptozotocin. The mechanism of action for streptozocin closely resembles that of 358 temozolomide as both react with DNA by undergoing substitution reactions forming a 359 methyldiazonium ion, resulting in methylation of primarily N 7 guanine (67%). They both 360 induce high levels of DNA methylation, and recognition and repair of this methylation 361 results in single-and double-strand DNA breaks. 29 To the best of our knowledge, no data 362 have been published on a correlation between hypermutation and streptozocin treatment, 363 but as streptozocin and temozolomide so closely resemble each other in their mechanism of 364 action, one can hypothesize the same mechanism to occur in streptozocin-treated patients. 365 It would be interesting to investigate whether prior treatment with streptozocin or 366 temozolomide indeed induces high TMB in mNEN, and if so, whether pretreatment with 367 streptozocin or temozolomide renders these tumors more sensitive to checkpoint inhibition. 368 Similarly, we observed a large contribution of the mutational signature associated with base 369 excision repair deficiency due to MUTYH aberrations in the second highest-TMB mNET, and 370 indeed this patient harbored a pathogenic germline MUTYH allele coupled with a complete 371 somatic loss of the respective chromosomal arm. MUTYH abnormalities have also previously 372 described to occur in pancreatic NET. 3 A single mNET presented a BRCA2-genotype 373 associated with homologous recombination deficiency but did not harbor (somatic) 374 mutations within BRCA2. It did harbor a somatic mutation in RAD51C, a gene known to be 375 involved with homologous recombination and repair of DNA. 376

377
Concerning genomic stability, we observed evidence of chromothripsis, a large-scale and 378 catastrophic chromosomal rearrangement, within six mNEN (four mNEC, two mNET). 379 Strikingly, four out of six chromothripsis events occurred on chromosome 12. In addition, 380 we observe the first occurrence of localized hypermutation (kataegis) in six mNEC. Kataegis 381 encompasses a pattern of localized hypermutations, which has been identified in various, 382 but not all and to a varying degree, cancer types. 30,31 These regions of kataegis often co-383 localize with regions of genetic rearrangements. Kataegis is thought to arise from frequent 384 genomic C-to-U deamination events as a result of APOBEC-family enzyme activity, a DNA 385 cytosine deaminase which was recently identified as an internal and thus far unrecognized 386 source of DNA damage and mutagenesis in various cancer types. 32 More recently, kataegis, 387 rather than tumor mutational burden, microsatellite instability or mismatch repair 388 P a g e 13 | 35 deficiency, was found to independently correlate with PD-L1/PD-L2 expression, and could 389 thus be a marker in response to immune checkpoint inhibition. 33 390 391 Using unbiased driver gene analysis (dN/dS and GISTIC2) on the mNEN cohort, and on 392 mNEC/mNEC separately to explore putative driver genes, we (re-)discovered 9 genes to be 393 structures the tumor can originate from (foregut, midgut or hindgut). When we compared 424 the various origins at a genomic level, we could observe an increasing TMB; ranging from 425 1.05 (mNET -Midgut; Q 1 -Q 3 : 0.75 -1.39) and 1.07 (mNET -Unknown; Q 1 -Q 3 : 0.84 -1.53) to 426 1.27 (mNET -Other; Q 1 -Q 3 : 1.10 -1.44) and 1.35 (mNET -Pancreas; Q 1 -Q 3 : 0.9 -2.12) to 5.45 427 (mNEC; Q 1 -Q 3 : 3.8 -8.85). In addition, when we compared the two largest groups of mNET 428 per primary localization (midgut and pancreas), we can readily distinguish between the two 429 subtypes based on somatic mutation and copy-number profiles. Yet strikingly, many midgut-430 derived mNET (n = 9; 23%) did not present a mutual driver gene but each was characterized 431 by distinct sets of mutated genes reflecting the heterogenous nature of the malignancy. were generated from 50-100 ng of genomic DNA using standard protocols (Illumina, San 482 Diego, CA, USA) and subsequently whole-genome sequenced in a HiSeq X Ten system using 483 the paired-end sequencing protocol (2x150bp) for both the metastatic tumor and matched 484 blood sample. 485 486 P a g e 17 | 35 Subsequent alignment, somatic mutation detection and in silico tumor cell percentage 487 estimation were performed in a uniform manner as detailed by Priestley et al. (2019). 14 488 Briefly, paired-end sequencing reads were aligned against the human reference genome 489 (GRCh37) using BWA-mem (v0.7.5a). 57 Duplicate reads were marked and small insertion and 490 deletions (InDels) were realigned using GATK IndelRealigner (v3.4.46). Prior to somatic SNV 491 and InDel variant calling, base qualities were recalibrated using GATK BQSR (v3.4.46). 58 492 Somatic SNV, InDels and MNV were called by Strelka (v1.0.14) using the matched peripheral representing the tumor purity-corrected variant frequencies, we followed a previously 529 described approach by Stephens et al. (2012) 68 , implemented as: 530 in which where is the ratio of primary-aligned and non-duplicated reads observed for 532 alternative allele m over the reference allele (VAF), is the in silico estimated tumor purity 533 fraction, is the absolute copy-number of the segment overlapping m and ℎ is the wild-534 type (healthy) copy number; ℎ = 2 for autosomes and allosomes in female samples and 535 ℎ = 1 for allosomes in male samples. 536

Detection and annotation of recurrent copy-number alterations 547
To detect recurrent copy-number alterations, we performed a GISTIC2 71 (v2.0.23) analysis 548 over the entire mNEN cohort and, again, four separate GISTIC2 analysis on the major 549 subgroups (mNEC, mNET and pancreas-and midgut-derived mNET). 550

551
The GISTIC2 was performed using the following settings: 552 Genes were annotated to GISTIC2 peaks (q ≤ 0.1) based on the following strategy; 553 1) GISTIC2 focal peaks (all_lesions.conf_95.txt) were overlapped to genes (from verified 554 and manually annotated loci, no pseudogenes or read-throughs and from standard 555 chromosomes; n = 36574) from GENCODE (GRCh37; v33), taking into consideration 556 only the genes overlapping with at least 100 base pairs within the detected GISTIC2 557 peak. 558 2) If a GISTIC2 focal peak overlapped with multiple GENCODE genes, a combined 559 database containing known drivers detected in a metastatic pan-cancer dataset 560  (supplementary figure 7d). By comparing the cophenetic 580 correlation coefficient, residual sum of squares and silhouette, we opted to generate seven 581 custom de novo signatures. Custom signatures were correlated to existing (COSMIC v3) 582 mutational signatures using cosine similarity. 583

Detection of chromothripsis 584
Shatterseek 18 (v0.4) using default parameters was used to detect chromothripsis-like 585 events. As input, we used the rounded absolute copy numbers (as derived by PURPLE) and 586 structural variants with an TAF ≥ 0.1 at either end of the breakpoint. The male sex 587 chromosome (chrY) was excluded. The criteria for a chromothripsis-like event were based 588 on the following criteria: a) total number of intra-chromosomal structural variants involved 589 in the event ≥25; b) max. number of oscillating CN segments (2 states) ≥7 or max. number of 590 oscillating CN segments (3 states) ≥14; c) total size of chromothripsis event ≥20 megabase 591 pairs (Mbp); d) satisfying the test of equal distribution of SV types (p > 0.05); and e) 592 satisfying the test of non-random SV distribution within the cluster region or chromosome 593 (p ≤ 0.05). 594

Classification of homologous recombination deficiency genotypes 595
To determine Homologous Recombination Deficiency (HRD) due to possible loss-of-function 596 of BRCA1 and/or BRCA2 (amongst others), we utilized the Classifier for Homologous 597 Recombination Deficiency with default settings (CHORD; v2.0). CHORD uses a random-forest 598 approach to classify samples into HR-deficient / HR-proficient categories. 78 599

Inventory of clinically-actionable somatic alterations and putative therapeutic targets 600
Current clinical relevance of somatic alterations in relation to putative treatment options or 601 resistance mechanisms and trial eligibility was determined based upon the following 602

Code availability 616
Analysis and visualization have been performed using the statistical platform language R 617 (3.6.2), all utilized custom code and scripts can be freely requested and distributed by 618 contacting the authors. 619  WGS biopsy (114x) and blood control (38x) n = 99 Multiple biopsies from same patient or incomplete records n = 13 Final inclusion of distinct metastatic biopsies (NEN) n = 86