Genetic heterogeneity of primary lesion and metastasis in small intestine neuroendocrine tumors

Data on intratumoral heterogeneity of small intestine neuroendocrine tumors (SI-NETs) and related liver metastasis are limited. The aim of this study was to characterize genetic heterogeneity of 5 patients with SI-NETs. Therefore, formalin-fixed, paraffin-embedded tissue samples of primary and metastatic lesions as well as benign liver of five patients with synchronously metastasized, well differentiated SI-NETs were analyzed with whole exome sequencing. For one patient, chip based 850k whole DNA methylome analysis was performed of primary and metastatic tumor tissue as well as control tissue. Thereby, 156 single nucleotide variants (SNVs) in 150 genes were identified and amount of mutations per sample ranged from 9–34 (mean 22). The degree of common (0–94%) and private mutations per sample was strongly varying (6–100%). In all patients, copy number variations (CNV) were found and the degree of intratumoral heterogeneity of CNVs corresponded to SNV analysis. DNA methylation analysis of a patient without common SNVs revealed a large overlap of common methylated CpG sites. In conclusion, SI-NET primary and metastatic lesions show a highly varying degree of intratumoral heterogeneity. Driver events might not be detectable with exome analysis only, and further comprehensive studies including whole genome and epigenetic analyses are warranted.

in patients with disease progression during first-line somatotstatin analogue therapy 9 . Development of targeted therapy approaches is impeded by only limited available data on genetic alterations and lack of potential driver genes of SI-NETs.
To date, three studies investigated the genomic landscape of SI-NETs with large scale next generation sequencing [10][11][12] . Thereby, SI-NETs were found to be genetically silent neoplasms with a mutational load of less than 1 mutation per mega base (Mb). Compared to the majority of other malignancies where mutation rates vary between 1 and >10/Mb, SI-NET can be regarded as one of the most silent neoplasms 13 . Data on variation of the mutational landscape between primary and metastatic lesions are limited and available exome data are highly diverse with the percentage of common mutations of primary lesion and metastasis ranging from 0-47% 11 . Further characterization of this intratumoral heterogeneity is of high importance to better understand tumorigenesis of SI-NETs and its possible implications on future therapy approaches.
In the present study we characterized the mutational landscape of pairs of primary tumor and hepatic metastases of well differentiated, synchronously metastasized SI-NETs from five patients. In addition, we performed DNA-methylome analysis in a patient with marked heterogeneity. We thereby unveiled a profound genetic variability strongly suggesting presence of subclonality within SI-NETs.

Results
Patients. Five patients with a SI-NET and synchronously diagnosed hepatic metastases were included. All tumors were positive for Synaptophysin, Chromogranin A, and CD56 and well/moderately differentiated (G1/ G2). Ki-67 was between <2% and 15. Estimated mean percentage of tumor cells was comparable for primary tumors (mean 80%, range 70-90%) and metastatic lesions (86%, 70-95%). One patient had a resection in curative intention but experienced recurrent disease after one year. At study closure, one patient had controlled tumor burden under therapy with somatostatin, in the other four patients 2 nd and 3 rd line therapies were applied due to progressive disease. Mean follow-up was 46 mo (25-65 mo) and 2/5 patients had passed away at study closure. A summary of clinical data is provided in Table 1.

Mutational analysis.
Exomes of all 15 tumor and non-tumor samples were successfully sequenced. A mean of 102 Mio reads (range 73-124 Mio) per sample was mapped. After removal of duplicates (mean 39%, range 31-54%), mean coverage depth was 120 reads per sample (range 87-146, Supplementary Table 2). Overall, 156 mutations in 150 genes were identified. Mean mutation rate per megabase was 0.36, while 0.31 mutations per megabase were identified in primary lesions and 0.40 in metastases, respectively, and amount of mutations per sample ranged from 9-34 (mean 22). 25 mutations were successfully validated with pyro-or Sanger sequencing (validation rate 100%). Functional categories and mutational characteristics are shown in Fig. 1A/B.

Somatic intratumoral heterogeneity.
We observed a markedly differing extent of intratumoral heterogeneity of primary tumor (p) and metastasis (m) in our cohort of SI-NETs. For example, in Patient 1 (P1) no common mutations were identified (0%). On the other hand, the amount of common mutations per sample varied among the other four patients such as 20% in P2m and 94% in P4p. Amount of private mutations only present in the primary lesion but not in the metastasis varied as well strongly in Patient 2-5 between 6% (P4p) and 57% (P2p). Likewise, private mutations only detected in metastasis but not in the primary tumor varied between 11% (P4m) and 80% (P2m). An overview of mutational distribution including representative histopathological images is provided in Fig. 2.
Besides mutational analysis, distribution of observed allele frequency (AF) was examined. Across all samples with enough variants (except the metastasis sample of patient 4 that did not harbour informative variants), we observed a markedly higher AF in metastasis samples than in the primary tumor. These differences were significant in 2/4 and a trend (p > 0.10) was observed in 1/4 tested patients (false discovery rate corrected t test; p < 0.020 for patient 1, p < 0.049 for patient 2, p < 0.104 for patient 4, p < 0.085 for patient 5, Fig. 1C).

Potential driver mutations.
To assess presence of potential driver mutations, SNVs in known cancer-associated genes or genes previously described to be mutated in SI-NET were classified as likely pathogenic or variants of unknown significance according to their SIFT score (see methods for details). Thereby, sixteen likely pathogenic SNVs and two variants of unknown significance were identified. Notably, in samples with overlapping SNVs (Patient 2-5), potential driver mutations were mainly found in both samples whereas Patient 1, who had no overlap of SNVs, harbored different likely pathogenic mutations in metastasis and primary lesion ( Table 2). The only recurrent non-synonymously mutated gene in this study was RBMS3, which was found being mutated in two patients. No marked difference in mutation pattern was observed: The most frequent SNV was C > T/G > A transition in both primary and metastatic lesions (Fig. 1B).
To investigate whether mutated functional clusters differ in primary lesion and metastasis, we compared mutated genes of primary SI-NET (including data of 103 primary lesions from previous studies) 10,11 and metastatic lesions (including data of five published metastatic lesions) 11 . We thereby found different canonical pathways to be mutated: For example, cell cycle proteins were identified in the primary lesions and gene<s of regulation of adherence junction stability and disassembly in metastasis (Supplementary Table 3).
Copy number variation analysis. The exome based copy number analysis revealed gains of chromosomes 4, 5, 7, 10, 14, 20 and 21 as well as losses of chromosome 13 and 18. Copy numbers of metastasis and primary were similar in Patient 3 and Patient 4, differed partially in Patient 2 and 5 and were completely different in Patient 1 (Fig. 3). Thereby, overlap of copy number variations (CNV) matched the overlap of SNVs in most patients (Patient 1, 2, 4, 5) but differed in Patient 3, where similar CNV alterations were observed despite heterogeneity in SNV analysis.

Methylation analysis.
To further investigate potential background of the differing mutational profiles of primary lesion and metastasis of Patient 1, additional analysis of whole DNA methylation status of both lesions was performed and compared to benign small intestine. Whole DNA-methylome data of average beta-values revealed markedly more epigenetic dysregulation in the metastasis compared to the primary lesion (Fig. 4A). Of interest, the majority of dysregulated CpG sites of the primary were found in the metastatic lesion as well: In Fig. 4B, overlap of all genes with differentially methylated, promotor associated CpG sites is shown. Within this overlap, one hypermethylated tumor suppressor gene was identified (KLF6). All differentially methylated genes are shown in Supplementary Table 4.
To investigate whether these alterations might be common within SI-NETs in terms of an entity specific signature, a comparison with a control cohort of methylation data of SI-NETs (n = 20) and benign small intestine (n = 20) was performed 12 . Thereby, the common methylation alterations of Patient 1 were not found in the control cohort (Supplementary Table 5).

Discussion
Characterization of the genetic landscape of SI-NETs is of high importance to identify potential targets for personalized therapy regimens. In the present study, primary and metastasis samples of five SI-NETs were investigated with whole exome sequencing to examine presence and extent of subclonality within this rare entity. Thereby, a profound genetic heterogeneity between primary lesions and hepatic metastasis was revealed.
A comparable mean somatic mutation rate of 0. 33  the comparability of our data with the mentioned studies, although DNA was extracted from formalin-fixed paraffin embedded (FFPE) tissue in the current study and not from frozen tissue as in the previous analyses.
While no relevant difference of the mutational pattern was observed between metastasis and primary lesion, mutational landscapes within our cohort differed clearly: A markedly varying amount of common mutations (present in both samples) as well as private mutations (not present in the corresponding sample) was identified. Moreover, in 1/5 cases no common mutations were observed at all. These observations correspond with data of Francis et al., where 2/3 of SI-NETs presenting with liver metastasis were found to have no common mutations in primary and metastasis 11 . In addition, AF-analysis revealed a markedly higher AF in metastasis samples than in the primary tumor, which could be explained by the metastasis growing out of a single clone from the primary while the primary has a more polyclonal population. These findings are highly remarkable since they indicate a different biological way of metastatic spread in comparison to other solid tumors. For example, recent studies on malignancies such as non-small cell lung cancer (>50% common mutations), endometrial cancer (>45%), and head and neck squamous cell carcinoma (>60%) observed high rates of common mutations in tumor and metastasis samples [14][15][16] . Moreover, in small cell lung cancer (SCLC), which is thought to originate from neuroendocrine cells as well, a rate of more than 95% common mutations was found suggesting a linear model of clonal evolution in this entity 14 . The observed heterogeneity has to be kept in mind for interpretation of previous genomic studies, which included only one sample of either primary or metastatic lesion. Furthermore, these data have to be considered in the planning of future biomarker-driven targeted therapy trials since results of a single biopsy might be insufficient. Copy number alterations are common events in SI-NETs. We also observed alterations in our cohort. For example, loss of chromosome 18 was found in 3/5 patients which was described to occur in >60% of SI-NETs. Likewise we found gains of chromosomes 4, 5, 7, 10, 14 and 20 as well as a loss of chromosome 13 [17][18][19][20] .
In-silico pathway analysis revelead that cell cycle pathways were found to be significantly mutated in the primary lesion, whereas genes regulating adherence junction stability and disassembly were identified in the metastases group. This correlates to protein expression data of a study of Kim et al., where proteins involved in cell cycle and proliferation were found to be differentially expressed in primary lesions and metastasis 21 . However, functional data of metastasis development of SI-NET are still due and further studies are warranted.
A possible explanation for the presence of genetically independent primary lesion and metastasis might be common epigenetic alterations, which could potentially lead to genetic instability and consecutive development of genetic subclonality. In line with this hypothesis, a recent study profiled 49 primary SI-NET tumors with a methylation array and identified a SI-NET specific panel of epigenetically altered genes with an associated change in the expression profile of the respective gene 12 . Even so, distinct epigenetic changes were observed as well in other tumors with a high degree of common mutations such as SCLC and a definite link between epigenetic changes and the SI-NET-specific heterogeneous mutational landscape has not been determined yet 22 . To investigate a potential common epigenetic alteration of the patient without common SNVs and CNVs (P1), we  Table 2. Potential driver genes. Identified potential driver mutations within the study cohort. Genes found to be mutated within this study present as well in the cancer gene census or previous studies on small intestine neuroendocrine tumors were analyzed for potential pathogenicity (see methods for details). Mutations with a SIFT score ≤0.05 as well as nonsense or frameshift mutations were categorized as likely pathogenic, whereas the rest was classified as variants of unknown significance (VUS). Pat: patient, NA: not available. performed an analysis of 850,000 methylation sites. We thereby found that the metastasis was clearly more epigenetically dysregulated in comparison to benign tissue than the primary lesion. which is in line with observations of a recent study investigating methylation profiles of SI-NET metastasis 23 . Moreover, we observed the majority of differentially (de)methylated promoter-associated CpGs of the primary to be present in the metastasis as well. Since these alterations were not found to be present in a control cohort, we hypothesize that these two lesions most likely originate from a common single clone. This suggests a common driver mutation, which is not detectable by exome sequencing. Since the majority of available sequencing data on SI-NETs is exome-based, these observations are of high interest. The only potential tumor suppressor gene being hypermethylated was KLF6. KLF6 depletion was associated with various cancers and a potential role in SI-NET warrants further investigation [24][25][26][27] . Unfortunately, since only FFPE tissue was available, no analysis of KLF6-expression could be performed within this study. Besides differential methylation, other epigenetic alterations such as histone modifications or mutations in regulatory elements have to be considered as potential underlying tumor drivers warranting further investigation. In addition, since epidemiological data suggest a higher risk of SI-NET development for first-degree relatives of patients with SI-NETs, germ-line alterations might predispose for tumor development 28 . This is strengthened by the observation that SI-NETs occur multiply in up to 30% 29 . Furthermore, a familial aggregation of SI-NETs with other noncarcinoid malignancies was reported 30 . However, to date, data of genome wide association studies on potentially predisposing polymorphisms for development of SI-NETs are lacking and further investigation is needed.
In conclusion, our data demonstrate a highly varying heterogeneity of SI-NETs. Moreover, tumor drivers of SI-NET development might not be detectable by exome sequencing and further comprehensive studies including whole genome sequencing as well as epigenetic analyses are needed.

Methods
Patients and sample preparation. Tissue and tumor samples as well as patient data used in this study were provided by the University Cancer Center Frankfurt (UCT). Written informed consent was obtained from all patients and the study was approved by the institutional Review Boards of the UCT and the Ethics Committee at the University Hospital Frankfurt (project-number: SGI-OW-01/2013). All methods were performed in accordance with the declaration of Helsinki. Only patients with available FFPE tissue of each the primary lesion, the metastasis and the matched normal tissue were included. The definite diagnosis of SI-NET was confirmed

Whole exome sequencing. Sequencing libraries were prepared from tumor and non-tumor tissue with
SureSelectXT Human All Exon V6 (target size 60 Mb, Agilent, Santa Clara, CA) according to the manufacturer's instructions and paired-end sequencing was performed on a HiSeq2500 (Illumina, San Diego, CA) with 2 × 100 base pairs (bp) read length. Input DNA was at least 1000 ng and at least 16 gigabases (range 16-20 Gb) of raw read data per sample were produced.

Variant calling.
Demultiplexing was performed with Illumina CASAVA (1.8.2) and adapters were trimmed with Skewer (0.1.116). Trimmed raw reads were aligned to the human genome (hg19) with the Burrows-Wheeler Aligner (BWA-mem version 0.7.2). Reads that aligned at more than one locus were discarded. Duplicate reads were removed with SAMtools (0.1.18). To enhance sensitivity, two different software tools were used for variant calling: SAMtools and varscan (2.3.5). Somatic variants in primary tumor and metastasis were selected by filtering against the variants present in the control sample. Then, all variants with a dbSNP ID were removed and only variants affecting coding exons were further evaluated. All called variants were manually reviewed in the Integrative Genomics Viewer (2.3) and all variants suspected of being technical artifacts were discarded. Only variants with an AF of ≥5% and at least ten reads were considered as true. For each mutation found only in one tumor sample, sequencing data of both the corresponding tumor and non-tumor sample were investigated for presence of reads with the same information. Alterations occurring only in either primary or metastasis were defined as private, whereas mutations present in both tumor sites were regarded as common. For analysis of AF, all variants with at least five reads were included. Pyrosequencing. For SNVs with low allele frequencies or in case primer construction for Sanger sequencing failed, pyrosequencing of both tumor samples and corresponding non-tumor sample was performed. Primer design was performed with PSQ Assay Design (Biotage, Uppsala, Sweden) and assays were established with Pyromark Q24 (Qiagen). PCR reaction was performed with the PyroMark PCR kit (Qiagen) according to manufacturer's recommendations using 20 pmol primer and 25-50 ng template DNA. The PCRs were performed as described for Sanger sequencing. The resulting PCR products were sequenced with the PyroMark Q24 pyrosequencer using PyroMark Gold Q96 reagents (Qiagen). All assays were run with tumor and non-tumor samples as well as positive (Qiagen human control DNA) and negative control. SNVs with ≥5% difference in mutant AF compared to non-tumor tissue and positive control were assessed as true variants.

Sanger sequencing. A representative selection of SNV and insertions and deletions (Indels
Identification of potential driver genes. All non-synonymous mutations detected in the current study which were present in the COSMIC cancer gene census (data extracted January 2017) were categorized as potential driver genes. Moreover, mutated genes were compared to former large scale sequencing studies on SI-NET, and recurrent genes were manually reviewed as well for potential oncogenic function 10,11 . Besides matching all non-synonymous mutations with the established list, all mutated genes were manually reviewed and those with high probability to have an oncogenic effect were examined in more detail. Mutations having passed this criteria with a SIFT score of ≤0.05, nonsense or frameshift mutations were classified as likely pathogenic mutations 31 . All other non-synonymous mutations of this list were classified as variants of unknown significance.
Copy Number analysis. CNV were computed on uniquely mapping, non-duplicate, high quality reads using an internally-developed method based on sequencing coverage depth. Briefly, we used at least 10 reference samples to create a model of the expected coverage that represents biases introduced by the target enrichment method, sequence GC content, library preparation protocol, insert size and sequencing technology, as well as inter-sample variation.
CNV calling for germline samples was performed by computing the sample's coverage profile, correcting for total read count and computing the deviation from the expected coverage. Genomic regions were called as variant if they deviated by at least 2 standard deviations from the model mean and the deviation was concordant with a biologically possible copy number (e.g., +50% for a heterozygous duplication, −50% for a heterozygous deletion). For tumor samples, the estimated tumor content was taken into account to deduct the copy number. For instance, given a tumor content of 60%, an observed deviation of +30% represented a heterozygous duplication in the tumor, while an observed deviation of +20% could either represent a heterozygous duplication of non-tumor tissue or a subclonal duplication in the tumor. To improve visual clarity and highlight large-scale changes, data was smoothed using the median over windows of five mega bases. Signal transduction pathway analysis. Functional cluster analysis was performed with the Genomatix Genome Analyzer (v3.70808, Genomatix GmbH, Munich, Germany) using the tool GeneRanker, which is based on FuncAssociate 32 . For pathway analysis, gene associations with over 400 canonical signal transduction pathways were performed collected from the Pathway Interaction Database and pathway commons 33,34 . All canonical pathways are derived from Homo sapiens. Genes were grouped into non-synonymously mutated genes in primary lesions as well as in metastasis of SI-NETs.

Methylation analysis.
Representative areas of tumor tissue were punched out of paraffin blocks. DNA was extracted using the Invisorb Genomic DNA Kit (Stratek, Birkenfeld, Germany). Amount of DNA was measured using Invitrogen's dsDNA BR Assay Kit and Qubit System (Invitrogen, Waltham, MA). Bi-sulfite treatment was performed using EZ DNA methylation Kit from Zymo Research (Zymo Research, Irvine, CA). DNA restoration was performed using DNA Clean and Concentrator-5 (Zymo Research) and Infinium HD FFPE Restore Kit (Illumina). All following steps were performed according to Infinium HD FFPE Methylation Assay manual protocol (Illumina). Pretreated DNA was processed and hybridized on an EPIC 850k chip (Illumina). Chips were scanned on an Iscan (Illumina). Data was processed using Illumina Genome Studio as well as JMP 11 (SAS). CpG-sites with potential SNPs were removed before analyses. First we compared average beta-values, ranging between 0 (fully unmethylated) and 1 (fully methylated) of primary tumor tissue, normal tissue, and a corresponding liver metastasis. Delta beta-values of whole methylome and promoter associated regions (as annotated by the manufacturer including transcriptional start sites such as TSS200, TSS1500) were compared for primary tumor vs. normal tissue and liver metastasis vs. normal tissue.
Methylation data were compared to a control cohort of already published methylation data of SI-NET primary lesions (n = 20) as well benign small intestine (n = 20) 12 .
Statistics. Descriptive statistics were calculated using BiAS (version 11.01, BiAS for Windows; Epsilon-Verlag, Frankfurt, Germany) and R 35 . P-values less than 0.05 were considered significant. Data availability. Exome data was submitted to the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI) with the accession number SRP126752. All other data are available on request from the corresponding author.