Longitudinal Copy-Number Alteration Analysis in Plasma Cell-Free DNA of Neuroendocrine Neoplasms is a Novel Specific Biomarker for Diagnosis, Prognosis, and Follow-up

Neuroendocrine neoplasms (NEN) are a heterogeneous group of neoplasms that develop from neuroendocrine cells present in various organs throughout the body, including pancreas, small intestine, and lung. NENs are historically considered rare with a combined age-adjusted incidence of 6.98 per 100,000 per year according to the Surveillance Epidemiology and End Results (SEER) database, with the highest incidence rates observed in gastroenteropancreatic NENs (GEP-NEN; 3.56 per 100,000) and lung NENs (1.49 per 100,000; ref. 1). However, their incidence has been increasing over the past decades and due to the frequently long disease course of NENs, GEP-NENs are now the second most prevalent malignancy of the digestive tract, after colorectal cancer (1). Currently, the only curative treatment option for NENs remains surgery, while only 50% of cases present with localized disease at time of diagnosis (1). The prognosis for patients with NEN is highly variable with median overall survival (OS) rates ranging from more than 30 years to less than 1 year and differs significantly according to primary site of the NEN (1). The World Health Organization (WHO) NEN classification system distinguishes well-differentiated NENs, i.e., neuroendocrine tumors (NET) that can be WHO grade (G) 1, 2, or 3 based on the proliferation rate, and poorly differentiated NENs, i.e., neuroendocrine carcinomas (NEC) that are more aggressive and always have a high proliferation rate (G3; ref. 2). Due to a low incidence and high clinical and biological heterogeneity, diagnosis of NENs remains challenging and depends mostly on (functional) imaging and histologic examination of tissue specimens. For pancreatic malignancies, it is crucial to distinguish pancreatic NENs (PNEN) from the more common pancreatic adenocarcinomas (PAAD), as this has important implications for prognosis and therapeutic management. Currently, the two most used circulating biomarkers in NENs are chromogranin A (CgA) and neuron-specific enolase (NSE). Although CgA is being considered the most useful circulating biomarker for management of NENs, it has several drawbacks including limited sensitivity and specificity, lack of assay standardization, and limited value as follow-up biomarker (3). NSE has poor diagnostic sensitivity and specificity, and although it might have some value as prognostic biomarker, it is generally considered as a suboptimal biomarker in NENs (4–6).

For many cancer types, liquid biopsies are emerging as minimally invasive biomarkers and especially the analysis of circulating cell-free DNA (cfDNA) is gaining a lot of interest. In patients with cancer, part of the cfDNA is tumor-derived, i.e., the circulating tumor DNA (ctDNA), which can be utilized as a tool for (repeated) molecular characterization of a tumor without requiring an invasive tissue biopsy (7). Copy-number alterations (CNA), for example, can be confidently detected in cfDNA in other tumor types, with a good correlation of detected CNAs between liquid and tissue biopsies (8, 9). Currently, knowledge and research on liquid biopsies in NENs remains limited to a few studies, despite the urgent need for new, preferably blood-based, biomarkers and the identification of several molecular characteristics that could aid in clinical management, but whose use is limited due to tissue availability (6, 10–15). CNAs represent an important molecular alteration in NENs and could thus be interesting to evaluate in cfDNA (12, 16–18). In addition, longitudinal monitoring of patients with NEN through sequential ctDNA evaluation also remains unexplored. Interestingly, CNA detection allows estimation of the fraction of tumor-derived DNA in total cfDNA, i.e., tumor fraction, which is considered as a marker for tumor burden and could be evaluated over time (7, 19).

The aim of this study was to explore biomarker potential of cfDNA CNA analysis in a prospective cohort of metastatic, well-differentiated lung and GEP-NENs, building upon our proof-of-concept study (10).

Within NETwerk, a network of eight Belgian hospitals [AZ Monica, AZ Klina, AZ Nikolaas, AZ Rivierenland, AZ Voorkempen, Gasthuiszusters Antwerpen, Antwerp University Hospital (UZA), and Ziekenhuis Netwerk Antwerpen], patients with well-differentiated, metastatic GEP-NEN and lung NEN were prospectively included at initiation of everolimus treatment and closely followed with 2- to 3-monthly imaging and (bi)monthly plasma collections, until progression, termination of everolimus treatment, or reaching the end of the study after 24 months.

In addition, plasma samples of patients with NEN were prospectively collected at the UZA, in collaboration with the Biobank@UZA (20), aiming for a 6-monthly sampling. Patients with a well-differentiated, metastatic GEP-NEN or a NEN of the respiratory system (lung NEN) of which at least one plasma sample was available were included in this study, including 2 patients from our previous study (10). All available, prospectively collected plasma samples of these patients were analyzed. Plasma samples were collected between 2013 and 2020.

In addition, plasma samples from 100 healthy, nonidentifiable subjects, i.e., females from routine noninvasive prenatal testing (NIPT), were used as reference samples.

Clinicopathologic data were collected from all patients, including age, sex, WHO 2017 NEN-grade, primary tumor site, tumor marker measurements, and survival data. CgA and NSE measurements performed within 28 days from plasma sampling were collected. OS was calculated as time since inclusion in our study (i.e., first sampling) until death or censoring. Progression-free survival (PFS) was calculated in everolimus-treated patients starting from initiation of everolimus treatment until investigator-assessed progression or censoring. In addition, "sum of longest diameters of target lesions" (DSUM) was determined according to RECIST 1.1.

The study was approved by the local ethics committees of all participating centers with the ethics committee of UZA/University of Antwerp as central ethics committee (approval numbers 17/28/316 and 16/46/490) and conducted in accordance with the Declaration of Helsinki. All participating patients provided written informed consent.

Publicly available whole-genome sequencing (WGS) data of 98 PNEN and 98 paired normal tissue samples (Dataset ID: EGAD00001002684; ref. 17) and WGS data of PAAD (N = 109) and normal reference (N = 40) tissue samples (Dataset ID: EGAD00001003927) were obtained. Both datasets were generated within the context of the International Cancer Genome Consortium and could be downloaded from the European Genome-Phenome Archive (EGA) after receiving approval from the dedicated Data Access Compliance Office (21).

In addition, we obtained shallow WGS (sWGS) data of 96 cfDNA samples of patients with PAAD, which were deposited at the NCBI Sequence Read Archive (SRA) under accession number PRJNA633741 (19).

Blood samples of everolimus-treated patients were collected in Cell-Free DNA BCT tubes (Streck, Biomedical Diagnostics, Antwerp, Belgium) and plasma was obtained using following two-step centrifugation procedure: (i) 1,600 g for 10 minutes and (ii) 16,000 g for 10 minutes. Plasma samples were then frozen at −80°C until further processing.

Blood samples of the biobanking project were collected in EDTA tubes and plasma was obtained after a centrifugation step of 10 minutes at 400 g or 1,500 g and frozen at −80°C. In addition, samples were centrifugated for 10 minutes at 16,000 g right before cfDNA extraction.

cfDNA was extracted from plasma using automated extraction with the QIAsymphony DSP Circulating DNA Kit on the QIAsymphony SP (Qiagen, Hilden, Germany). The Qubit dsDNA High Sensitivity Assay Kit for the Qubit 2.0 fluorometer (Thermo Fisher Scientific) was used to measure cfDNA concentrations.

sWGS was performed on the cfDNA samples. Library preparation was performed using the TruSeq Nano DNA Library Prep Kit (Illumina) on an automated Hamilton STAR Liquid Handling System (Hamilton Germany GmbH – Robotics, Gräfelfing, Germany) with 2.5 to 12.5 ng of cfDNA as input. Then, single-end 75 bp sequencing of the samples was performed on the NextSeq 500/550 (Illumina). The whole workflow was optimized and validated in context of diagnostic NIPT for chromosomal abnormalities at the Center of Medical Genetics Antwerp.

Raw reads were mapped to human reference genome hg19 using the Burrows-Wheeler Aligner (BWA) minimum essential medium (MEM) algorithm (v0.7.4; ref. 22). Duplicate reads were removed by Picard MarkDuplicates (v2.21.7; http://broadinstitute.github.io/picard/). SAMtools (v1.9) was used for sorting and indexing of the BAM files (23). CNA calling and estimation of tumor fraction were performed using the R-based tool ichorCNA (8). The input files for ichorCNA are wig files with read count data. These were generated using the readCounter function of the HMMcopy Suite by counting the number of reads that were mapped per 50 kb nonoverlapping window or bin in which the genome was divided (24). Only reads with a mapping quality above 20 were withheld by readCounter. To improve our ichorCNA analysis, 100 control samples from nonidentifiable, healthy subjects were used to generate a “Panel of Normals” (PoN) reference. Additionally, the nontumor fraction start values were set to 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, and 0.99, and to increase the specificity of segmentation the txnE and txnStrength parameters were set to 0.9999999 and 10,000,000, respectively. Only autosomal chromosomes were analyzed.

The genome-wide guanine-cytosine (GC) and mappability corrected mean absolute deviation (MAD) value of copy ratios between adjacent bins was calculated by ichorCNA for every sample. According to the ichorCNA manual, MAD values below 0.150 and coverages higher than 0.1x indicate high-quality data, which we therefore used as cutoffs. Data interpretation was performed according to the ichorCNA manual and ichorCNA solutions were visually checked for intrapatient consistency. The tumor fraction was estimated for every sample collection timepoint of a patient and a tumor fraction greater than 0.03, the cutoff determined by ichorCNA, indicated the presence of ctDNA. To determine the tumor fraction evolution during follow-up, a representative CNA profile as determined by ichorCNA was selected for every patient and the tumor fraction for this specific pattern was selected at every timepoint.

The WGS data of the PNEN tissue samples (Dataset ID: EGAD00001002684) and PAAD tissue samples (Dataset ID: EGAD00001003927) were available as BAM files (17, 21). The tumor samples were first downsampled to approximately 5x using SAMtools and subsequently analyzed for CNAs using the same approach as the cfDNA samples, i.e., readCounter (of the HMMcopy Suite) followed by ichorCNA analysis (8, 23, 24). The normal, reference samples were used to generate a PoN for ichorCNA analysis, and as the tumoral fraction in tissue was generally higher than in cfDNA, nontumor fraction start values were 0.3, 0.5, and 0.8, and default specificity values were used. In accordance with cfDNA, only autosomal chromosomes were analyzed.

In addition, the sWGS data of 96 cfDNA samples of patients with PAAD (NCBI SRA: PRJNA633741) could be downloaded as raw fastq files using the fasterq-dump function of the SRA Toolkit v2.9.6 (https://github.com/ncbi/sra-tools; ref. 19). We have only used the first read of the paired-end sequencing, and these were analyzed with exactly the same approach and settings as our own cfDNA samples. As for this cohort, no reference samples analyzed with the same protocol were available; we used our own cfDNA PoN.

The ichorCNA output of all samples was checked for an MAD < 0.150 and a coverage > 0.1x, indicating high-quality samples.

A classification model was generated to predict the type of pancreatic malignancy, i.e., PNEN or PAAD, based on the observed CNA pattern. PNEN and PAAD tissue samples were randomly distributed in a training set (80%) and a validation set (20%) and their ichorCNA log_R values were used as input data. First, feature extraction was performed using principle component analysis (PCA) on the training set and the first 10 principle components (PC) were then used to generate a binary logistic regression classification model. Next, only the significant PCs were selected to generate the final binary logistic regression classification model. To apply the model on new data, PCs were calculated in the same way as those of the training set whereafter the model could be applied on the calculated PCs for classification. The generated model was used to predict the class of all tissue samples in the training and validation set. Additionally, the classification model was tested on cfDNA samples with detectable ctDNA by classifying them as originating from a patient with PAAD or a patient with PNEN. Sensitivity, specificity, and area under the receiver operating characteristic curve (AUC) for every classification were determined. The classification modeling was performed in R version 3.6 and the R packages “ROCR” and “InformationValue” were used (25).

All statistical analyses were performed in R version 3.6. Associations between two categorical variables were examined using Fisher exact test. The t test was used to compare the means of normally distributed continuous variables between two groups. CgA, NSE, and cfDNA concentration were first log-transformed to obtain a normally distributed variable, after which the t test could be applied. Correlation testing for nonnormal distributed continuous variables was performed using Spearman correlation analysis.

OS analysis was performed via univariate and multivariate Cox proportional hazards (CoxPH) regression modeling and proportionality was confirmed using Schoenfeld residuals. Due to lower number of events for OS analysis in the subset analysis per location of primary tumor, association was examined using the log-rank test. The survival curves were generated using the Kaplan–Meier method.

In a subpopulation of patients undergoing everolimus treatment and regular sample collections, an association between tumor fraction evolution and PFS was evaluated using a joint modeling approach. A joint model combines longitudinal biomarker measurements and time-to-event data in a single model to evaluate an association. The R package “JM” was used to perform the joint modeling (26). A joint model was fitted for the data using a random slope model with only time as fixed effect and a Cox survival model with single intercept. The method used to generate a joint model determines the type of survival submodel to be fitted and the type of numerical integration, and we have used the piecewise-PH-aGH method which uses the relative risk model with a piecewise-constant baseline.

All P values were based on two-sided hypothesis testing and the applied cutoff for statistical significance was 0.05.

The data can be obtained from the corresponding author upon reasonable request.

In total, 195 longitudinal plasma samples (median, 3 samples/patient; range, 1–21 samples/patient) were prospectively collected from 43 patients, including both patients with GEP-NEN (N = 40) and lung NEN (N = 3; Table 1). Two patients had an insulinoma, while all other patients had a nonfunctional NEN. Two patients were diagnosed with the MEN1 syndrome. Since there were no inclusion restrictions regarding treatment, patients received a variety of treatments during sample collection (Table 1).

cfDNA was extracted from all 195 plasma samples, subjected to sWGS, and analyzed using ichorCNA. One sample did not meet the quality criteria and was removed for further analyses. The median coverage, calculated after deduplication and quality filtering, was 0.3x (range, 0.2x–0.4x). The analysis detected ctDNA in 58 of the 194 cfDNA samples (30%) with a median tumor fraction of 0.15 (range, 0.031–0.92). In 19 out of 43 patients (44%), at least one sample was ctDNA positive (ctDNA⁺).

All 100 control samples from healthy individuals were ctDNA negative (ctDNA⁻), validating the predetermined tumor fraction cutoff for our cohort.

Associations between ctDNA positivity and clinicopathologic characteristics were examined (Table 1). ctDNA positivity was significantly associated with higher WHO grade (P = 0.037), presence of primary tumor in situ (P = 0.034), and location of the primary tumor (P = 0.011). A relatively high proportion of PNENs (62%) and both the stomach and primary tumor unknown (PTU) NEN were ctDNA⁺, while only 24% of small intestinal NENs (siNEN) were ctDNA⁺ and all 3 lung NENs were ctDNA⁻. Treatment at inclusion was significantly associated with ctDNA positivity (P = 0.018), with patients receiving a later line of treatment, such as everolimus, being more frequently ctDNA⁺.

The mean cfDNA concentrations (in ng/mL plasma) were significantly different between ctDNA⁻ and ctDNA⁺ samples, with higher cfDNA concentrations measured in ctDNA⁺ samples (P = 6.6e-09; Supplementary Fig. S1A). Furthermore, mean CgA values (N = 92) were significantly different between ctDNA⁺ and ctDNA⁻ samples, with higher values observed for ctDNA⁺ samples (P = 0.00059; Supplementary Fig. S1B). In total, 51% of samples were concordant for CgA and ctDNA, e.g., both positive or both negative, 46% were CgA⁺/ctDNA⁻, while 3% were CgA⁻/ctDNA⁺. The CgA⁻/ctDNA⁺ samples originated from 2 patients for whom ctDNA analysis could thus increase the diagnostic yield compared with CgA analysis alone. In addition, mean NSE values (N = 55) were significantly different between ctDNA⁺ and ctDNA⁻ samples, with higher values observed for ctDNA⁺ samples (P = 0.0076; Supplementary Fig. S1C). In total, 55% of the samples were concordant for NSE and ctDNA, 40% were NSE⁺/ctDNA⁻, and 5% were NSE⁻/ctDNA⁺. The only sample of 1 patient and two samples of another patient were NSE⁻/ctDNA⁺, ctDNA analysis might thus provide additional information for these patients/timepoints.

In our cohort, median follow-up time was 27 months (range, 0.39–89.7 months). The median time between inclusion and the first ctDNA⁺ sample was 0 months (range, 0–42.8 months), with 12 patients already having a ctDNA⁺ sample at inclusion. Univariate CoxPH modeling for ctDNA positivity showed that ctDNA⁺ patients had a significantly higher risk of death compared with ctDNA⁻ patients, with an HR of 5.6 [95% confidence interval (CI): 1.5–20, P = 0.009; Fig. 1A]. In a multivariate analysis including ctDNA status and location, encoded as pancreatic versus extrapancreatic, ctDNA positivity was significantly associated with a higher risk of death with an HR of 4.5 (95% CI, 1.17–17.3, P = 0.029), while tumor location was not significant (HR: 1.9; 95% CI, 0.57–6.7, P = 0.29).

Since location of the primary tumor was significantly associated with ctDNA positivity and is also known to determine prognosis of the patients with NEN, we performed a subset analysis based on tumor location. In the PNEN subgroup, the association between ctDNA positivity and OS remained significant (log-rank P = 0.0089; Fig. 1B). Median OS in the ctDNA⁺ group was 17 months (95% CI, 10.7–not reached) and was not reached (95% CI, not reached–not reached) in the ctDNA⁻ group. An association could not be observed in the subset of siNENs (log-rank P = 0.45, median OS not reached in both groups).

The CNA profile in tumoral tissue of 98 patients with PNEN was determined, showing detectable CNAs in all samples and a median tumor fraction of 0.729 (range, 0.154–0.999; ref. 17). Frequency plots were generated for the CNAs detected in cfDNA of our ctDNA⁺ patients with PNEN and the PNEN tissue dataset (Fig. 2). Upon visual inspection of the frequency plots, the CNA patterns in cfDNA and tumoral tissue were very comparable and in addition, they were also significantly correlated (Spearman correlation coefficient = 0.82, P < 2.2e-16). This result validates the PNEN tumoral origin of the detected CNA patterns in cfDNA.

Since the detected CNAs have a specific pattern in PNENs, we hypothesized that the detected CNA pattern could distinguish PNENs from the more common PAADs. Classification of tumor types based on CNA patterns has been previously done for various tumor types, but not yet for PNENs and PAADs (27). To investigate this, a classification model was first generated and validated on a tissue dataset of PNENs and PAADs, and afterwards tested in a PNEN and PAAD cfDNA dataset. Therefore, CNA profiles of PAAD tissue and cfDNA samples were determined, and a tissue sample without detectable CNAs was removed for further analyses. Frequency plots for CNAs in PAAD tissue (N = 108) and ctDNA⁺ cfDNA samples (N = 84) are shown in Supplementary Fig. S2.

The final binary logistic regression classification model, consisting of five PCs, had a sensitivity of 94% and specificity of 95% with an AUC of 99% for the tissue training set. Sensitivity, specificity, and AUC in the validation cohort were 100%, 95%, and 100%, respectively.

The model was then used to predict the tumor type of all ctDNA⁺ PNEN (N = 45) and PAAD (N = 84) cfDNA samples. In the cfDNA cohort, the classification model reached a sensitivity of 62% and specificity of 86% with an AUC of 79%. The quality of the classification slightly decreased when performed on cfDNA samples compared with tissue with a loss in sensitivity, but still a good specificity. Both samples of patient P05 were classified as PAAD samples, while a tumor biopsy of this patient taken at inclusion showed dedifferentiation of a G3 NET towards a G3 NEC in the biopsied lesion. In addition, samples from patient P07, who had a 5-HIAA–producing PNEN, were also consistently misclassified.

In our study, 18 patients (P01–P18) were prospectively included at initiation of everolimus treatment after disease progression and subjected to a standardized follow-up. Upon visual inspection of the CNA patterns of ctDNA⁺ samples within patients, little changes were observed over time. In patient P01, the only observed difference between the first and the last sample was a subclonal gain in chromosome 15 (Fig. 3A). In P04, several additional subclonal gains were observed at one of the timepoints; however, the tumor fraction was very high at this timepoint, which leads to a higher sensitivity for subclonal alterations and therefore does not necessarily reflect newly developed alterations. In P05, more pronounced changes were observed between the two samples, which were taken 3 weeks apart (Fig. 3B). This could be explained by the dedifferentiation and progression under treatment that was observed in this patient. In P06, a subclonal loss of the long arm of chromosome 6 had emerged 4 months after inclusion.

The date of disease progression after everolimus treatment initiation was determined for all patients and PFS was calculated. To evaluate a possible association between longitudinal measurements of tumor fraction and PFS, a joint model was fitted. The generated model had a significant association value of 3.6137 (P = 0.035), which is a measure for the effect of the tumor fraction in the risk for progression. An increase in tumor fraction resulted in a decreased PFS probability (Fig. 4).

Figure 5A and B illustrate the tumor fraction evolution over time in 2 representative patients. In P06, after reaching the lowest tumor fraction since inclusion, the tumor fraction increased for 4 months after which progressive disease was also observed on imaging. In P17, the tumor fraction dropped to undetectable levels after treatment initiation and durable stable disease was observed on imaging in this patient.

In Fig. 5C, the tumor fraction evolution during long-term follow-up of P04 is shown. The first ctDNA⁺ sample was observed during a treatment break due to selective coronary artery bypass grafting surgery, but the tumor fraction afterwards decreased to undetectable levels. Later, the tumor fraction started to steeply increase, which overlapped with an increase in tumor size (DSUM) and finally also investigator-assessed disease progression. A switch to CAPTEM treatment resulted in a drastic drop in tumor fraction, followed by a decrease in tumor size. The decrease in tumor fraction might thus be an early indicator of response to treatment.

The first longitudinal study evaluating cfDNA in blood in a cohort of well-differentiated, metastatic lung and GEP-NENs is presented here. Biomarker potential for cfDNA regarding diagnosis, prognosis, and follow-up was explored by generating a genome-wide overview of tumor-associated CNAs detectable in cfDNA of 195 prospectively collected plasma samples from 43 patients with NEN that were compared with healthy control cfDNA (N = 100). All plasma samples derived from patients with PNEN (N = 21) were used for comparison with publicly available PNEN tissue (N = 98), PAAD tissue (N = 109), and PAAD cfDNA (N = 96) data.

To this day, the main origin of cfDNA remains uncertain as several biological processes are likely to contribute to the release of cfDNA in the bloodstream (28). It has been suggested that cfDNA enters the bloodstream through two different processes, namely passive and active release. The first process refers to cell death via both apoptosis and necrosis, while the second process involves the release of extracellular vesicles containing DNA. Generally, patients with cancer present with elevated levels of cfDNA, which could be explained by an increase in cell death in large, advanced tumors (29). This has been partially confirmed in our study since cfDNA levels were higher in ctDNA⁺ samples. In total, 30% of the analyzed NEN plasma samples and 44% of the patients had detectable ctDNA. Highly variable ctDNA detection rates have been observed in different types of advanced cancers ranging from 10% to 100% of cases, which includes our detection rate (30). A potential explanation for this variability could be that growth kinetics differ between tumor types. For example, colorectal cancer has a cell loss factor of 96% and generally presents with higher levels of ctDNA. This suggests that high cell losses due to cell death could potentially lead to increasing release of ctDNA in the blood stream (31). To our knowledge, these kinetic parameters have not yet been determined in NENs, but it is known that NENs are generally slow growing and therefore may have a lower cell loss factor, hence relatively low ctDNA release. This could explain why ctDNA was not detected in all patients. This hypothesis is strengthened by the fact that the sample with the highest ctDNA concentration originated from a patient who had just received a cell death–inducing ablation which may have caused this increase in the amount of ctDNA. All our 100 analyzed control samples from healthy individuals were ctDNA⁻, making this analysis highly specific for ctDNA detection.

First, we compared clinicopathologic characteristics between ctDNA⁺ and ctDNA⁻ patient groups. ctDNA positivity was significantly associated with higher WHO grades, which is in line with findings for other tumor types. In lung cancer, for example, it has been described that tumors with a higher proliferation rate also have a higher rate of ctDNA shedding (32). Location of the primary tumor was also significantly associated with ctDNA positivity with a large proportion of PNENs being ctDNA⁺ (62%), while only 24% of siNENs were ctDNA⁺ and all 3 lung NENs were ctDNA⁻. These differences in ctDNA positivity might be explained by differences in tumor biology between NENs of different primary origin, which could also result in differences in ctDNA shedding. In addition, ctDNA⁺ samples had a significantly higher cfDNA concentration and the CgA and NSE values were also significantly higher around ctDNA⁺ timepoints. Although sensitivity might be lower than currently used tumor markers, such as CgA and NSE, ctDNA analysis is highly specific, which is very important for rare diseases. The lower specificity is an important drawback of the current tumor markers, e.g., a CgA increase can be caused by several other conditions, including renal and hepatic dysfunction, cardiovascular problems, and use of proton pump inhibitors, leading to false-positive results (3). Additionally, ctDNA analysis provides additional information in selected patients and timepoints with normal CgA and NSE values.

The so-called “NETest”, developed by Modlin and colleagues, is a new, circulating biomarker assay for NENs based on the detection of a panel of circulating transcripts in the blood. Quantitative analysis of these transcripts allows the calculation of a disease-activity score, which could be useful for diagnosis, follow-up, and prediction of Peptide Receptor Radionuclide Therapy (PRRT) treatment response (33). However, the results of additional, independent validation studies and practical aspects, including availability, cost, and transparency, will further determine the integration of the NETest in clinical practice. An important advantage of our cfDNA analysis over the NETest is the potential to obtain a global profile of molecular alterations in the tumor, compared with a numeric score for the NETest. In addition to gathering genome-wide information, the technique used in our article has several other advantages including high-throughput potential, cost-effectiveness, standardization, availability, and easy implementation as many labs are already performing it in context of NIPT analysis. Efforts will be made to continuously improve our strategy. The cost of our analysis is already much lower as compared with the standard follow-up scans, and we will try to further decrease this cost during further optimizations. The same applies to the high-throughput potential. Currently, several samples are run simultaneously. However, the number of samples as well as the speed of the analysis will be improved in order to accelerate the turnaround time, which is currently 1 week. Our findings should also be confirmed in a larger cohort of patients with NEN as this analysis was performed in a rather limited cohort. Another drawback of our current cfDNA analysis is the limited sensitivity observed in our cohort. Recently, Zakka and colleagues have performed targeted sequencing of cfDNA in a retrospective and heterogeneous NEN cohort in which they illustrated the feasibility of ctDNA evaluation in NENs, detecting alterations in 87.5% of patients (11). Several aspects could contribute to the difference in ctDNA detection with our study. First, targeted sequencing has a higher sensitivity compared with CNA detection with sWGS, allowing detection of alterations up to fractions of 0.0002. However, this sensitivity also allows detection of variants linked to "clonal hematopoiesis with indeterminate potential" that can be misinterpreted as tumor-related variants, such as the commonly implicated gene TP53, and for which they did not correct (11, 34, 35). In addition, up to 48% of their patients with NEN were NECs, which have higher proliferation rates and are expected to shed more ctDNA (32). Additional drawbacks were the use of a general gene panel, which lacked genes with an important role in NENs such as DAXX/ATRX/MEN1 in PNENs, and the very limited clinical data, which is crucial for data interpretation and translation towards the clinic (17, 36). The sensitivity for CNA detection of our cfDNA analysis might thus be increased by a paired, deeper, or targeted sequencing approach, which improves sensitivity and resolution and allows detection of alterations on gene level (37). A promising targeted sequencing approach that has been successfully applied for CNA detection is based on the use of single-molecule molecular inversion probes (smMIP) and next-generation sequencing. This strategy was already used in studies on sensitive and reliable molecular diagnostics for patients with cancer (38, 39). Consequently, this seems an extremely interesting approach to potentially improve the sensitivity of our CNA analysis by combining sWGS with the use of smMIPs. Increased sensitivity might also be achieved by combining CNA detection with mutation and/or epigenetic profiling of cfDNA, and could be interesting to explore further.

With this prospective study, we demonstrate that tumor-associated CNA detection is feasible in longitudinal plasma cfDNA samples. This represents a first step in the development of improved liquid biomarkers for identifying, prognosticating, and monitoring progression in patients with NEN, an important unmet need according to the European Neuroendocrine Tumor Society (ENETS; ref. 40). These liquid biomarkers could be extremely useful in patients where a tissue biopsy cannot be performed or histologic examination is inconclusive. CNA patterns detected in cfDNA of our patients with PNEN and in public PNEN tissue were very comparable, validating the tumoral origin of the CNAs we detected in cfDNA (17). In most patients, the CNA pattern in cfDNA also remained stable during our follow-up period. It is of course still possible that during the complete disease course, which can be very long for these patients, changes do occur.

Interestingly, classification of tumors based on their tissue CNA patterns has been performed for other tumor types, as CNA patterns were found to differ between different tumor types (27). However, besides a study in patients with lung cancer, classification based on plasma cfDNA CNA patterns remains unexplored (9). Since we detected quite specific CNA patterns in PNENs, we explored the possibility of a classification model that could distinguish PNENs from the more common PAADs using CNA patterns. The model reached AUCs of 99% and 100% in the tissue training and validation set, respectively. Classification of 84 PAAD cfDNA samples and our 45 PNEN cfDNA samples resulted in a sensitivity, specificity, and AUC of 62%, 86%, and 79%, respectively (19). Despite the slightly decreased performance of the classification model for cfDNA samples, it already demonstrates diagnostic potential for classification of PNENs versus PAADs using cfDNA CNA patterns with good specificity. As observed, classification might be affected by cases with more uncommon molecular alterations that are underrepresented in the training set, like 5-HIAA–producing PNENs, or by NET dedifferentiation. NET dedifferentiation could be associated with changes in CNA patterns, but molecular knowledge regarding this process remains scarce so far (41). Additional research is thus required to further explore the diagnostic potential of cfDNA analysis and should focus on extension of the training and validation cohorts and the use of more advanced classification techniques, such as neural networks, for generation of the most optimal classification model.

Univariate OS analysis in the whole cohort indicated a significantly worse OS in the group of ctDNA⁺ patients. In addition, the difference remained significant in multivariate analysis including location and subset analysis in patients with PNEN. Patients that become ctDNA⁺ during their disease course thus have a worse prognosis, and cfDNA analysis might therefore be useful as a prognostic biomarker. This prognostic potential is promising especially in patients with PTU NENs. The sole patient with PTU NEN included in our analysis appeared to be ctDNA⁺ which implied, based on the OS analysis, that this patient had a significantly worse OS as compared with ctDNA⁻ NENs. This is in line with previous findings which revealed that patients with PTU NENs have poor prognoses (42). Our findings should be confirmed in larger cohorts of patients with NEN.

Everolimus is a frequently used targeted therapy for patients with NEN, but in line with other targeted therapies, primary and acquired resistance to everolimus poses an important limitation for its application. Unfortunately, we are only beginning to understand molecular resistance mechanisms and predictive markers are still lacking (43, 44). Therefore, timely detection of resistance is crucial and is currently mainly based on imaging, which is not always straightforward. Therefore, we have explored tumor fraction measurements for follow-up of everolimus-treated patients. In our study, 18 patients with NEN were included at everolimus treatment initiation and followed via a standardized protocol, including frequent plasma sampling and imaging. A joint modeling approach showed a significant association between longitudinal tumor fraction measurements and the risk for progression. The model predicted a decreased PFS probability with increasing tumor fraction levels. A joint model could ultimately lead to individual survival predictions based on tumor fraction measurements, but further research in a more extensive cohort is therefore required. Furthermore, we have illustrated tumor fraction evolution in relation to the patients' disease course. Hereby, increasing tumor fractions could be observed before disease progression, while a decreasing tumor fraction could be observed in patients with durable stable disease. Patient follow-up using tumor fraction measurements might thus be interesting to evaluate response or progression under treatment, and could assist in making treatment decisions.

In conclusion, we demonstrated biomarker potential for longitudinal plasma cfDNA CNA analysis in a cohort of lung and GEP-NENs with multiple clinical applications. We illustrated that the detected CNAs in cfDNA were tumor-derived and that CNA patterns in cfDNA could assist in diagnostic classification of PNENs and the more common PAADs. Furthermore, patients with detectable ctDNA had a worse OS and longitudinal tumor fraction measurements were associated with PFS and might be useful in anticipating tumoral progression. Our research is thus an important first step towards clinical implementation of ctDNA analysis for NENs and shows promise for diagnosis, prognosis, and follow-up and as an alternative to solid biopsies for detection of tumor-related molecular alterations.

G. Boons reports grants from Kom op tegen Kanker and Research Foundation - Flanders during the conduct of the study. T. Vandamme reports grants from Kom op tegen Kanker during the conduct of the study. L. Mariën reports grants from Kom op tegen Kanker during the conduct of the study. M. Peeters reports grants from Kom op tegen Kanker during the conduct of the study. No disclosures were reported by the other authors.

G. Boons: Conceptualization, resources, formal analysis, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. T. Vandamme: Conceptualization, resources, formal analysis, investigation, visualization, methodology, project administration, writing–review and editing. L. Mariën: Validation, writing–review and editing. W. Lybaert: Resources, investigation, writing–review and editing. G. Roeyen: Resources, investigation, writing–review and editing. T. Rondou: Resources, investigation, writing–review and editing. K. Papadimitriou: Resources, investigation, writing–review and editing. K. Janssens: Resources, investigation, writing–review and editing. B. Op de Beeck: Resources, investigation, writing–review and editing. M. Simoens: Resources, investigation, writing–review and editing. W. Demey: Resources, investigation, writing–review and editing. I. Dero: Resources, investigation, writing–review and editing. G. Van Camp: Conceptualization, resources, investigation, methodology, project administration, writing–review and editing. M. Peeters: Conceptualization, resources, investigation, methodology, project administration, writing–review and editing. K. Op de Beeck: Conceptualization, resources, investigation, methodology, project administration, writing–review and editing.

This work was supported by Kom op tegen Kanker (Stand up to Cancer, the Flemish Cancer Society). G. Boons is supported by a Ph.D. fellowship of the Research Foundation – Flanders (FWO; 1195118N). We would like to thank all participating centers and investigators for their contribution to this research. In addition, we would like to thank Lesley De Backer and Isolde Van der Massen from the Multidisciplinary Oncology Center of Antwerp (MOCA) and the NETwerk for their support in collecting clinicopathological data, Ella Roelant for her support in statistical analysis, and our bioinformatics team for their support. Further, we would like to thank the lab technicians of the NIPT team for their contribution.

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.