Comparative Targeted Genome Profiling between Solid and Liquid Biopsies in Gastroenteropancreatic Neuroendocrine Neoplasms: A Proof-of-Concept Pilot Study

Gastroenteropancreatic neuroendocrine tumors (GEP-NETs) originate from the diffuse neuroendocrine system of the gastrointestinal tract and pancreas. They represent 1–4% of all gastrointestinal neoplasms, and their incidence is recording a considerable worldwide increase due to advances in diagnostic techniques and increased awareness, with a European prevalence ranging from 2.1 to 6.6 cases per 100,000 population [1, 2]. They present a highly heterogeneous clinical behavior, varying according to site of origin, differentiation grade, proliferative rate, and functioning properties [3]. A patient-tailored management is desirable to differentiate those patients who could benefit from more conservative treatment from those who need early aggressive therapy. Therefore, predictive and prognostic tools are needed. Literature data agree about the prognostic, diagnostic, and predictive role of Ki-67 as an expression of the tumor proliferation rate determining the tumor grade (G1–G2–G3) [4, 5]. Other useful characteristics that should be considered in the prognostic evaluation are tumor size, disease extension, differentiation, functional status, tumor staging, patient age, SSTR immunostaining, and molecular imaging [4, 6, 7]. Clinical scores were also found helpful in stratifying survival probability [8, 9]. GEP-NET genetic profile could represent another tool to better characterize tumor behavior and outcome, and recent improvements in next-generation sequencing (NGS) technologies have enhanced the exploration of GEP-NET genetic background [10, 11]. Scarpa et al. [11] described the mutational landscape of pancreatic NETs (pNETs), based on the whole-exon sequencing of 102 primary pNETs. They found common somatic mutations in genes involved in chromatin remodeling and telomere maintenance (such as DAXX, ATRXX, SETD2, ARID1, and MEN1), DNA damage repair (such as MUTYH and CHECK2), activation of mTOR signaling (such as TSC1, TSC2, MTOR, PTEN, DEPDC5) [10, 11]. Over the last few years, liquid biopsy (LB) for the detection and analysis of circulating tumor DNA (ctDNA) has become an attractive alternative to traditional tissue solid biopsies (SBs) for diagnostic, prognostic, therapeutic, and monitoring purposes [12] resulting more cost effective, less invasive, and potentially more representative of tumor heterogeneity [13, 14]. To date, ctDNA analysis from LB has rarely been applied to NETs since recurrent mutations have not been frequently identified, differently from other tumors [12, 15, 16]. Furthermore, comparative analysis focusing on SB and LB comparison is lacking. Therefore, the aim of our pilot study was to compare the mutational profile (MP) provided by LB with that deriving from SB in GEP-NETs.

This multicenter study enrolled patients with a diagnosis of GEP-NETs who were addressed to primary tumor/metastasis surgery as standard treatment practice in two Italian centers (Azienda Ospedaliero, Universitaria of Ferrara, and IRCCS Ospedale Policlinico San Martino of Genova). Blood samples were collected at the time of surgery. Each patient who accepted to participate in the study provided written informed consent. The Declaration of Helsinki guidelines and the International Conference on Harmonization (ICH) Good Clinical Practice (GCP) guidelines were followed for ethical conduct during the study.

Both liquid and tissue samples underwent high-throughput targeted-exome sequencing to evaluate the presence of mutations in 11 genes previously reported as altered in NETs [10, 11]. Concerning LB, peripheral blood samples were collected and processed within 2 h from withdrawal for each patient. After plasma isolation, samples were stored at −80°C until DNA extraction. The QIAamp Circulating Nucleic Acid Kit and the QIAamp DNA FFPE Tissue Kit from Qiagen (Hilden, Germany) were employed following the manufacturer’s protocol to extract DNA from LB and SB, respectively. Prehybridization library preparation was performed using the Sureselect low-input LP kit from Agilent Technologies (CA, USA) with minor modifications. Up to 10 μL of each sample was taken as input and adjusted to a volume of 50 μL using molecular-grade water. Enzymatic fragmentation was omitted, and the samples were directly subjected to end repair by incubating at 20°C for 15 min and then at 72°C for 15 min. Following end repair, unique molecular identifiers were added. The samples were purified using SPRI beads and quantified using the Qubit 2.0 BR Assay from Thermo Fisher (MA, USA). Additionally, all samples underwent MiSeq Nano (Illumina, CA, USA) analysis to quantify the prehybridization generated libraries before equimolar pooling. A maximum of 500 ng/sample was combined in a single pool. Hybridization was carried out using custom idt probe ID 04022561 and the xGen Hybridization and Wash Kit from Integrated DNA Technologies (Coralville, IA, USA). The final libraries were quantified using the Qubit 2.0 DNA HS assay and TapeStation HS D1000 Screen Tape from Agilent Technologies (CA, USA). The library quantity was measured using KAPA SYBR FAST qPCR with the QuantStudio 5 System from Applied Biosystems (California, USA). Equimolar pooling of libraries was performed based on qPCR QC values prior to sequencing on an Illumina NovaSeq S4 (Illumina, CA, USA) with a read length configuration of 150 paired-end for 13.3 million paired-end reads (6.75 million in each direction) per sample. Targeted-exome sequencing of selected genes (Table 1) was conducted on the Illumina NovaSeq S4 following the manufacturer’s protocol. Genetic analyses have been performed by an external laboratory (Biodiversa s.r.l.; Treviso, Italy).

The FASTQ files were processed to obtain paired reads with unique molecular identifiers using the “UMI_tool” and its “extract” function. After renaming the FASTQ files, the paired reads were trimmed using “TRIMMOMATIC” with quality and length parameters set at 20 and 35, respectively. The trimmed reads were then mapped to the hg19 reference genome. Subsequently, the resulting BAM files were sorted, and duplicate reads were marked using “sambamba markdup.” BAM tags (XS = sample, XC = cell, RX = umi) were assigned to the resulting BAM files using the bamtag utility from the “umis” tool. Reads originating from the same original molecule were grouped using the GroupReadsByUmi utility. Consensus sequences were identified using the CallMolecularConsensusReads utility and filtered using the FilterConsensusReads utility, both requiring a minimum of 3 reads to form and filter consensuses. The GroupReadsByUmi, CallMolecularConsensusReads, and FilterConsensusReads are utilities from the “fgbio” package. The resulting consensus reads were then remapped using bwa-mem to generate analysis-ready BAM files. Read sequences were extracted using the software “samtools fastq” to perform variant calling analysis, which was performed using the Sention^® Genomics software. Paired reads were mapped to the reference genome, and TNhaplotyer was used for somatic variant calling analysis of the tumor, using a reference panel of normal data to remove germinal variants. Annotation of variants was conducted using ANNOVAR, a pipeline for comprehensive functional annotation of genomic variants. OMIM information was added to the VCF file, and variants were filtered based on coverage (DP >6). Only variants that passed the quality control (filter PASS) were considered eligible for further analysis.

VCF files were converted into annovar VCF format using the wANNOVAR web server (https://wannovar.wglab.org/#get-started). All downstream analyses were performed in R version 4.2.1 within the RStudio environment. The annovar VCF files were then transformed into Mutation Annotation Format files using the R package maftools and merged into a single comprehensive Mutation Annotation Format file. Subsequent analyses were conducted on the entire dataset as well as for the comparison between LB and SB. For each group, the following analyses were performed: classification of variants, type of variants, class of single-nucleotide variants (SNVs), variants per sample, summary of variant classification, percentage of variants per genes for each group, ratio of transitions to transversions in single-nucleotide polymorphisms, ratio of transition to transversion analysis for individual genes, and variant allele frequency analysis.

Categorical data were summarized using frequencies and percentages and were compared using Fisher’s exact test. Continuous variables were reported as median and were compared by the Mann-Whitney U test. A p value <0.05 was considered significant. Hierarchical clustering was performed separately for each gene using a data matrix to define a more specific dataset for hierarchical clustering analysis. Euclidean distance-based hierarchical clustering was conducted for each gene. The distance matrix was constructed based on the number of exonic mutations per sample.

Six patients were included in the study, two women and four men, with a median age of 64.5 years (range 54–74). At enrollment, they were all affected by a metastatic nonfunctioning GEP-NET; in particular, 3 patients were diagnosed with an ileal NET and 3 patients with a pNET. One patient harbored a MEN1 germline mutation. Baseline characteristics of the whole cohort are summarized in Table 2. For each patient, we collected tissue samples from the primary tumor (Ge03p, Ge04p, Ge06p) or metastatic site (Fe02p, Fe05p, Fe06p), and we obtained a corresponding LB (50_s, 43_s, 44_s, 45_s, 46_s, 41_s) at the time of surgery. Plasma samples were collected the day before surgery.

A comparison analysis was performed between SB and LB, and then, for each patient, a head-to-head comparison was performed. The analysis revealed 262 mutations in SB and 271 mutations in LB, with a median number of variants/sample of 55 and 66.5, respectively. SB and LB were similar as concerns number and type of mutations (Table 3). Conversely, more significant differences resulted after performing a sample-to-sample comparison (Fig. 1). Total number of mutations was strictly sample-dependent for both SB and LB, as well as the LB_{(total mutations)}/SB_{(total mutations)} ratio. SB and LB prevalent mutations were represented by missense mutations and nonsense mutations, mainly determined by single-nucleotide polymorphisms. The highest percentage of variants found consisted of a C>T transition (Fig. 2). Concerning the mutational burden, SB presented 0.2–1.9 mut/Mb and LB 0.3–1.8 mut/Mb. In both cases, the translational effect was mostly non-synonymous (Fig. 3). Among the investigated genes, only ARID1A was mutated in all patients considering both SB and LB. All solid samples also harbored MTOR and ATRX gene alterations. On the other hand, all LB displayed MEN1, TSC2, PTEN, SETD2, and MUTYH gene alterations. CHEK2 was less frequently mutated, and DAXX mutations were not found in LB (Table 4). To explore the concordance between the MP in LB and SB samples, the latter were compared within each patient (Table 5). This analysis revealed 17 recurrent mutations that were shared between SB and LB, except for the patient FE06p/44_s, who has no shared mutations between SB and LB samples (Table 6). Interestingly, this patient is the only one bearing a germline MEN1 mutation. Two SNVs of the MTOR gene (c.G4731A and c.C2997T) were shared by 5 out of 6 patients, while MTOR SNV c.T1437C and SETD2 SNV c.T3333C were shared by 4 out of 6 patients. MEN1 SNV c.T1314C, MTOR SNV c.T1437C, and PTEN SNV c.G194C were shared among 2 out of 3 primary tumor samples; furthermore, two silent MTOR mutations, c.G4731A and c.C2997T, were also detected in 2 out of 3 metastatic samples, alongside the MTOR SNV c.T1437C. Finally, a hierarchical clustering analysis was performed on a customized data matrix built using exonic non-silent genetic mutations, and clinical data from each sample. The cluster dendrogram generated from the binary matrix revealed distinct clustering patterns among the samples. Notably, the Ge06p/41_s and Fe05p/43_s samples cluster together using this approach, indicating a high level of similarity between these samples in terms of their MP. On the other hand, the remaining samples displayed lower genomic similarity (Fig. 4).

Our study focused on evaluating a set of gene mutations previously found in GEP-NETs [10, 11], comparing the MP assessed in SB and in corresponding LB from GEP-NET patients, exploring the feasibility of LB evaluation for tumor characterization purposes. The interest in LB is increasing because it represents an attractive alternative to traditional tissue biopsies for diagnostic, prognostic, therapeutic, and monitoring purposes. Indeed, LB is easy to perform with low patient risks, could be easily repeated during follow-up, and could be used to evaluate treatment response. In addition, LB potentially reflects tumor heterogeneity, especially in metastatic disease [12‒14], allowing for ctDNA isolation for molecular profiling. ctDNA represents only 0.1–10% of the total circulating cell-free DNA, and its levels might depend on tumor load, stage, and response to treatment. Recently, a multigenomic blood-based assay has been proposed, the NETest [15‒18]. This tool may be useful for diagnostic [15, 19], prognostic [20, 21], and predictive purposes in NET patients [21‒23]. However, van Treijen et al. [24] showed a decrease in NETest score prognostic accuracy over time. Retrospective studies comparing mutation detection in plasma against that in matched tumor samples reported a sensitivity for ctDNA analysis between 65 and 98% [12]. To date, analysis of ctDNA from LB has rarely been applied to GEP-NET, since they present a relative heterogeneity in recurrent mutations in comparison with other tumors, lacking a reliable “mutational signature” [12, 15]. Knappskog et al. [15] compared the SB and LB sequencing of 76 cancer-related genes in 50 patients with advanced GEP neuroendocrine carcinomas, finding a concordance in 71% of cases. Furthermore, ctDNA concentration demonstrated the ability to discriminate metastatic and localized pNETs and to predict progression-free survival [12, 19]. LB and SB have been previously explored in a NET cohort including patients with medullary thyroid cancer or lung, pancreas, thymus, parathyroid, stomach, and rectum localizations reporting actionable mutations [25]. However, they mixed somatic and germline investigations as well as LB, SB, and blood samples, hampering a comparison with our study. Indeed, to our knowledge, our pilot study is the first to perform a concordance analysis between DNA sequencing from SB and LB samples in a cohort of GEP-NET patients. Each sample was tested for a panel of selected genetic mutations frequently associated with pNETs [10, 11] to achieve a high mutational rate that could help the comparative analysis between SB and LB sequencing. Our cohort included 3 stage IV pNETs in progression at evaluation, and SB was taken at metastatic sites. The other 3 patients were affected by a metastatic ileal NET, and SB was obtained from the primary tumor site at the time of surgery. Previous literature data have already shown high ctDNA levels in patients with pNETs and in ileal NETs compared to healthy controls [12, 19]. In addition, metastatic disease at evaluation could have been responsible for a high ctDNA concentration in our cohort, possibly influencing mutation detection. First, we compared MP derived from SB to that from LB, obtaining similar results in terms of number, type of mutations, and mutational burden. These findings suggest that LB MP is comparable to SB MP. On the other hand, C>T transitions were frequent in SB but not in LB, where mostly transversion was detected. The latter finding could possibly indicate a more heterogeneous and unstable disease [26, 27], reflecting genetic alterations developing during the metastatic process. In addition, in both SB and LB, the translational effect was mostly non-synonymous, suggesting a functional impact on corresponding proteins. Moreover, all evaluated genes were frequently mutated, suggesting their contribution in ileal NETs, besides their possible role as driver genes in pNETs. Only DAXX mutations were completely absent in LB. ARID1, DEPDC5, and CHECK2 mutation rates were the same in SB and LB. For other genes, LB mostly showed higher mutational rates, supporting ctDNA as a mirror of the complexity of the disease. Conversely, sample-to-sample comparison showed striking differences. Except for patient Fe02p/50_s, SB and LB deriving from the other patients showed a pronounced imbalance in mutational rates. In particular, LB 43_s, 45_s, and 41_s showed more mutations as compared to the SB counterparts (Fe05p, Ge03p, and Ge06p). The opposite occurred for SB and LB from patients Fe06p/44_s and Ge04p/46_s. These differences were associated neither with tumor and SB sites nor with tumor grade, stage, and proliferative activity. These results support the hypothesis of an important variability in NET genetic background that are reported to lack a strong mutational signature [12].

The concordance analysis between SB and LB MP demonstrated good performance, except for patient FE06p/44_s, who had no shared mutations. This finding may be due to the reported DNA damage accumulation and genome instability in conditionally MEN1 mutated mice [28]. Metastatic sites may shed a high ctDNA amount, possibly masking the detection of the mutations displayed by the seeding cancer cells. In addition, this patient was later on found to have a slowly progressing immunoproliferative disease, contributing in the peculiar heterogeneity of the mutational landscape in this patient. This hypothesis is further supported by the evidence that Ge04p/46_s and Fe05p/43_s samples cluster together at hierarchical analysis, indicating a high level of similarity between these samples in terms of their MP. Nine of 17 shared mutations involved genes contributing to the PI3K/AKT/mTOR pathway, supporting the key role of PI3K/AKT/mTOR signaling hyperactivation in NETs and the rationale of mTOR inhibitor therapeutic use. Finally, our results support the applicability of LB in GEP-NET molecular profiling.

In conclusion, our data show that tissue-based NGS and plasma-based NGS display a comparable sensitivity in the detection of clinically relevant mutations in GEP-NETs. Although preliminary, our findings strongly encourage the use of plasma-based multigene NGS panels as powerful and less invasive method to assess GEP-NET MP. In addition, LB results have potential therapeutic applications, especially in those settings where SB is difficult to obtain, such as lung neuroendocrine tumors. The latter may display actionable mutations (i.e., in receptor tyrosine kinases) that may be identified by LB, and the patient may be consequently addressed to a potentially effective treatment [29, 30].

Gastroenteropancreatic neuroendocrine neoplasm is a complex, heterogeneous, and multifactorial pathology resulting from the alteration of multiple regulatory pathways. In this frame, poor and heterogeneous sampling may negatively affect the stability and precision of the mutational scenario [31, 32], thus preventing the discovery of diagnostically functioning biomarkers. Therefore, the methodology and analytical approach employed in our study may arouse concern. Indeed, the statistical power of the present study is low, given the largely suboptimal sample size and the heterogenicity of the samples. Consequently, the mutational scenarios herein described, as well as the potential biomarkers, should not be considered as having diagnostic significance.

Notwithstanding the above-described statistical limitations, the main aim of the present pilot study, namely, exploring the mere possibility of using LB in GEP-NETs for MP evaluation, has to be considered fulfilled, while other studies, with high statistical accuracy and power, are needed to evaluate stability and precision. The main limitation of this study is indeed the small and heterogeneous sample size, which hampers the generalizability of our findings and prevents to draw broad conclusions that are beyond the scope of our findings. In addition to these limitations, a different tissue origin, disease stage, grade, and duration represent additional potential biases possibly influencing the interpretation of our results. However, as a pilot study using a consolidated multigenic panel targeted for pNET driver genes, it paves the way to future prospective and longitudinal studies that will analyze larger and homogeneous cohorts of different NET subtypes to evaluate LB clinical reliability and applicability at diagnosis, during follow-up, and treatment response evaluation. Larger studies with more homogeneous patient populations, focusing on specific NET subgroups. are indeed warranted to provide deeper insights into the differences and clinical applicability of LB as compared to SB.

All patients provided written informed consent to participate in the study. Ethics approval was not required due to the paucity of eligible patients.

The authors have no conflicts of interest to declare.

The study was funded by the PRIN grant 2017-PRN-PR.A-ZM_001 from MIUR and by FAR 2022 and 2023 from the University of Ferrara.

Conceptualization: M.C.Z., A.F., and M.A.; methodology: L.C.; investigation: I.G., P.B.S., L.R., F.C., M.M., R.M., A.L.S., and F.G.; formal analysis: L.C., I.G., P.B.S., and L.R.; writing – original draft preparation: I.G. and P.B.S.; writing – review and editing: M.C.Z., L.R., and M.M.; supervision: M.C.Z. All authors read and approved the final manuscript.

Irene Gagliardi and Federica Campolo contributed equally to this work.

The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.