Abstract
Purpose: Neuroblastoma displays important clinical and genetic heterogeneity, with emergence of new mutations at tumor progression.
Experimental Design: To study clonal evolution during treatment and follow-up, an innovative method based on circulating cell-free DNA (cfDNA) analysis by whole-exome sequencing (WES) paired with target sequencing was realized in sequential liquid biopsy samples of 19 neuroblastoma patients.
Results: WES of the primary tumor and cfDNA at diagnosis showed overlap of single-nucleotide variants (SNV) and copy number alterations, with 41% and 93% of all detected alterations common to the primary neuroblastoma and cfDNA. CfDNA WES at a second time point indicated a mean of 22 new SNVs for patients with progressive disease. Relapse-specific alterations included genes of the MAPK pathway and targeted the protein kinase A signaling pathway. Deep coverage target sequencing of intermediate time points during treatment and follow-up identified distinct subclones. For 17 seemingly relapse-specific SNVs detected by cfDNA WES at relapse but not tumor or cfDNA WES at diagnosis, deep coverage target sequencing detected these alterations in minor subclones, with relapse-emerging SNVs targeting genes of neuritogenesis and cell cycle. Furthermore a persisting, resistant clone with concomitant disappearance of other clones was identified by a mutation in the ubiquitin protein ligase HERC2.
Conclusions: Modelization of mutated allele fractions in cfDNA indicated distinct patterns of clonal evolution, with either a minor, treatment-resistant clone expanding to a major clone at relapse, or minor clones collaborating toward tumor progression. Identification of treatment-resistant clones will enable development of more efficient treatment strategies. Clin Cancer Res; 24(4); 939–49. ©2017 AACR.
In cancer with genetic heterogeneity, the identification of clones resistant to upfront treatment will lead to new therapeutic strategies targeting these resistant clones. We present a new approach based on whole-exome sequencing (WES) techniques of cell-free DNA (cfDNA) samples in neuroblastoma patients, enabling for the first time sequential analyses of mutational profiles in these patients. We now show how the study of cfDNA by WES can contribute to the understanding of clonal evolution in this malignancy, indicating the emergence of subclonal events present at diagnosis to clonal events at disease progression. The modelization of clonal evolution highlights mechanisms of treatment resistance and of escape. The identification of treatment-resistant clones will enable adaptation of treatment strategies.
Introduction
Analysis of circulating tumor DNA (ctDNA), a fraction of cell-free DNA (cfDNA), is a revolutionary tool for the study of tumor-specific genetic alterations in patients with cancer. Extracted from blood, it can be used as a surrogate marker for molecular diagnosis, and for estimation of tumor burden (1–3). These techniques enable an access of tumor molecular information when a biopsy is impossible and render sequential molecular analyses possible. Shift of mutational patterns over time has been demonstrated in particular in leukemia, but in most solid tumors, sequential analysis of tumor tissue is not possible due to the limited number of accessible tumor tissue time points (4–8).
Neuroblastoma, the most frequent extracranial solid tumor of early childhood, is characterized by both clinical and genetic heterogeneity with few recurrent biomarkers. Amplification of the MYCN oncogene is observed in 20% of cases and is associated with poor prognosis. Other copy number alterations (CNA) occur over more extensive chromosome regions, with the segmental chromosome alterations being associated with a poor outcome (9). Few recurrent mutations have been described, the most frequent targeting ALK (10%–12%), genes involved in chromatin remodeling (ATRX, ARID1A), or TERT rearrangements.
Both spatial and temporal genetic heterogeneity play an important role in neuroblastoma. Spatial heterogeneity has been described within a tumor, and between a tumor and its metastasis, for genetic features such as MYCN, or other CNAs (10). Temporal heterogeneity has also been observed with an accumulation of new CNAs, and mutations including ALK and RAS-MAPK mutations at the time of progression (11–13). However, more extensive studies of these phenomena are hampered by the scarcity of paired diagnosis–relapse samples, highlighting the importance of surrogate markers.
In neuroblastoma, several studies have indicated the presence of ctDNA in blood, enabling the detection of MYCN or ALK alterations (14–16). The overall copy number profile can also be determined in cfDNA using different technologies such as commercially available platforms or low coverage whole-genome sequencing (WGS; refs. 10, 17, 18).
The characterization of both CNA and a full mutational spectrum, of importance when analyzing mechanisms of tumor progression or treatment resistance, would be possible by whole-exome sequencing (WES) of cfDNA. Although a few limited studies have demonstrated its feasibility, the proposed protocols have not been applied to larger scale studies (19–21).
We now demonstrate for the first time the feasibility of WES cfDNA based on a modification of standard WES protocols enabling sequencing of extremely low amounts of cfDNA, enabling both CNA and SNV analysis in sequential neuroblastoma cfDNA samples. Deep coverage target sequencing was performed for validation and for analysis of intermediate time points (Fig. 1).
Materials and Methods
Patients and samples
Patients with neuroblastoma were included in this study if samples for WES of the primary tumor and plasma obtained at diagnosis were both available (n = 19 patients; Table 1; Supplementary Table S1). For 2 patients, plasma samples obtained at diagnosis were available, 9 had plasma collected at diagnosis and relapse, and 9 had 2 to 6 plasma collected during treatment and follow-up (Supplementary Fig. S1).
Case . | Clinical information . | Plasma: cfDNA Quantity (ng/mL) . | Genetic alterations detected at the time of diagnosis (0- none; 1-primary tumor only; 2-cfDNA only; 3-both primary tumor and cfDNA) . | ||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
. | INSS Stage at diagnosis . | Age at diagnosis (M) . | Delay diagnosis-2nd time point . | Status at 2nd time point . | Delay diagnosis-last follow-up (M) . | Status at last follow-up . | At diagnosis . | At 2nd time point . | Genomic copy number profile . | MYCN Amplification . | 1p- . | 1q+ . | 2p+ . | 3p- . | 4p- . | 11q- . | 14q- . | 17q+ . | Genes known to be recurrently altered in neuroblastoma or in other tumors . |
1 | 4 | 20 | 7 | SD | 14 | Dead | 13,533 | 71 | Numerical | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | MAP3K1: 3; MLLT3: 3 |
2a | 3 | 4 | 53 | CR | 47 | Alive | 406 | 104 | Segmental | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | MLLT4: 3; RUNX3: 2; ATXNL3: 2; GATA6: 2 |
3 | 4 | 25 | 7 | CR | 50 | Alive | 2,170 | 193 | Segmental | 3 | 3 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | ADCY4: 3; KRAS: 1; STAG2: 1 |
4 | 3 | 19 | 6 | CR | 51 | Alive | 76 | 262 | Segmental | 3 | 3 | 0 | 0 | 3 | 0 | 0 | 0 | 3 | ALK: 3; MAP2K2: 3 |
5 | 4 | 47 | 6 | CR | 40 | Alive | 1,295 | 1267 | Segmental | 0 | 0 | 0 | 0 | 3 | 3 | 3 | 0 | 3 | ARHGAP29: 3; ARHGAP33: 3 |
6 | 4 | 101 | 1 | PR | 7 | Dead | 794 | 371 | Segmental | 3 | 3 | 3 | 0 | 0 | 0 | 0 | 3 | 0 | ARHGAP18: 3; ITIH2: 3 |
7 | 4s | 3 | 7 | PD | 16 | Dead | 1,248 | 3617 | Segmental | 3 | 3 | 3 | 0 | 0 | 0 | 0 | 0 | 3 | |
8 | 4 | 19 | 12 | CR | 27 | Alive | 1,113 | 232 | Segmental | 3 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | MLL4: 2; NBPF10: 3 |
9 | 3 | 20 | 4 | CR | 36 | Alive | 2,328 | 73 | Segmental | 3 | 2 | 0 | 2 | 0 | 0 | 0 | 0 | 2 | ARHGAP6: 2; ERG: 2 |
10 | 4 | 32 | 27 | Relapse | 31 | Dead | 387 | 368 | Segmental | 0 | 0 | 0 | 3 | 3 | 0 | 3 | 0 | 3 | NKTR: 2; RASA1: 2; FLG:3 |
11 | 4 | 11 | 23 | Relapse | 31 | Dead | 639 | 105 | Segmental | 0 | 3 | 3 | 0 | 0 | 0 | 3 | 0 | 3 | ALK: 3; AXIN1: 3; MAP2K7: 3; NOTCH1: 3 |
12 | 1 | 58 | 8 | Relapse | 69 | Alive | 210 | 137 | Segmental | 0 | NC | ITPR2: 3 | |||||||
13 | 4 | 31 | 22 | Relapse | 24 | Dead | 1,278 | 364 | Segmental | 3 | 3 | 3 | 0 | 0 | 0 | 3 | 0 | 3 | ADCY4: 3; ADCY7: 3; ARHGAP30: 3; CHD5: 3; XRCC5: 3; RAD23A: 3 |
14 | 4 | 14 | 9 | Relapse | 11 | Dead | 930 | / | Segmental | 0 | 3 | 3 | 3 | 0 | 0 | 0 | 0 | 3 | ALK: 3; MAPK10: 1; BRCA2: 1; CDK19: 3 |
15 | 4 | 30 | 17 | Relapse | 18 | Dead | 1,530 | 376 | Segmental | 0 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | ARHGAP18: 2; DICER1: 2; NOTCH2:2; SHH: 2 |
16 | 4 | 66 | 26 | Relapse | 49 | Alive | 39 | 65 | Segmental | 0 | 3 | 3 | 2 | 3 | 3 | 3 | 0 | 3 | NOTCH4: 3 |
17 | 2 | 30 | 14 | Relapse | 21 | Dead | 26 | 30 | Segmental | 0 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | RNASEL: 3 |
18 | 2 | 9 | / | / | 12 | Dead | 350 | / | Segmental | 0 | 3 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | CDK3: 3 |
19 | 4 | 18 | 14 | Relapse | 21 | Dead | 42 | 1,210 | Segmental | 3 | 3 | 0 | 0 | 0 | 0 | 3 | 0 | 3 | MLL5: 1 |
Case . | Clinical information . | Plasma: cfDNA Quantity (ng/mL) . | Genetic alterations detected at the time of diagnosis (0- none; 1-primary tumor only; 2-cfDNA only; 3-both primary tumor and cfDNA) . | ||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
. | INSS Stage at diagnosis . | Age at diagnosis (M) . | Delay diagnosis-2nd time point . | Status at 2nd time point . | Delay diagnosis-last follow-up (M) . | Status at last follow-up . | At diagnosis . | At 2nd time point . | Genomic copy number profile . | MYCN Amplification . | 1p- . | 1q+ . | 2p+ . | 3p- . | 4p- . | 11q- . | 14q- . | 17q+ . | Genes known to be recurrently altered in neuroblastoma or in other tumors . |
1 | 4 | 20 | 7 | SD | 14 | Dead | 13,533 | 71 | Numerical | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | MAP3K1: 3; MLLT3: 3 |
2a | 3 | 4 | 53 | CR | 47 | Alive | 406 | 104 | Segmental | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | MLLT4: 3; RUNX3: 2; ATXNL3: 2; GATA6: 2 |
3 | 4 | 25 | 7 | CR | 50 | Alive | 2,170 | 193 | Segmental | 3 | 3 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | ADCY4: 3; KRAS: 1; STAG2: 1 |
4 | 3 | 19 | 6 | CR | 51 | Alive | 76 | 262 | Segmental | 3 | 3 | 0 | 0 | 3 | 0 | 0 | 0 | 3 | ALK: 3; MAP2K2: 3 |
5 | 4 | 47 | 6 | CR | 40 | Alive | 1,295 | 1267 | Segmental | 0 | 0 | 0 | 0 | 3 | 3 | 3 | 0 | 3 | ARHGAP29: 3; ARHGAP33: 3 |
6 | 4 | 101 | 1 | PR | 7 | Dead | 794 | 371 | Segmental | 3 | 3 | 3 | 0 | 0 | 0 | 0 | 3 | 0 | ARHGAP18: 3; ITIH2: 3 |
7 | 4s | 3 | 7 | PD | 16 | Dead | 1,248 | 3617 | Segmental | 3 | 3 | 3 | 0 | 0 | 0 | 0 | 0 | 3 | |
8 | 4 | 19 | 12 | CR | 27 | Alive | 1,113 | 232 | Segmental | 3 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | MLL4: 2; NBPF10: 3 |
9 | 3 | 20 | 4 | CR | 36 | Alive | 2,328 | 73 | Segmental | 3 | 2 | 0 | 2 | 0 | 0 | 0 | 0 | 2 | ARHGAP6: 2; ERG: 2 |
10 | 4 | 32 | 27 | Relapse | 31 | Dead | 387 | 368 | Segmental | 0 | 0 | 0 | 3 | 3 | 0 | 3 | 0 | 3 | NKTR: 2; RASA1: 2; FLG:3 |
11 | 4 | 11 | 23 | Relapse | 31 | Dead | 639 | 105 | Segmental | 0 | 3 | 3 | 0 | 0 | 0 | 3 | 0 | 3 | ALK: 3; AXIN1: 3; MAP2K7: 3; NOTCH1: 3 |
12 | 1 | 58 | 8 | Relapse | 69 | Alive | 210 | 137 | Segmental | 0 | NC | ITPR2: 3 | |||||||
13 | 4 | 31 | 22 | Relapse | 24 | Dead | 1,278 | 364 | Segmental | 3 | 3 | 3 | 0 | 0 | 0 | 3 | 0 | 3 | ADCY4: 3; ADCY7: 3; ARHGAP30: 3; CHD5: 3; XRCC5: 3; RAD23A: 3 |
14 | 4 | 14 | 9 | Relapse | 11 | Dead | 930 | / | Segmental | 0 | 3 | 3 | 3 | 0 | 0 | 0 | 0 | 3 | ALK: 3; MAPK10: 1; BRCA2: 1; CDK19: 3 |
15 | 4 | 30 | 17 | Relapse | 18 | Dead | 1,530 | 376 | Segmental | 0 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | ARHGAP18: 2; DICER1: 2; NOTCH2:2; SHH: 2 |
16 | 4 | 66 | 26 | Relapse | 49 | Alive | 39 | 65 | Segmental | 0 | 3 | 3 | 2 | 3 | 3 | 3 | 0 | 3 | NOTCH4: 3 |
17 | 2 | 30 | 14 | Relapse | 21 | Dead | 26 | 30 | Segmental | 0 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | RNASEL: 3 |
18 | 2 | 9 | / | / | 12 | Dead | 350 | / | Segmental | 0 | 3 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | CDK3: 3 |
19 | 4 | 18 | 14 | Relapse | 21 | Dead | 42 | 1,210 | Segmental | 3 | 3 | 0 | 0 | 0 | 0 | 3 | 0 | 3 | MLL5: 1 |
NOTE: cfDNA WES at diagnosis only permitted SNV calling but not CNA analysis.
Abbreviations: SD, stable disease; M, months; NC, noncontributive.
aCase 2: stage IVs at diagnosis with serum available at diagnosis. The study of this case is based on a local relapse with WES on tumor and plasma at relapse and target sequencing on plasma follow-up.
Patients were treated in French centers of the SFCE (Société Française de lutte contre les Cancers et les leucémies de l'Enfant et de l'adolescent) according to the relevant protocols (Supplementary Information S1). Written informed consent was obtained from parents/guardians according to national law. Samples were collected following inclusion of patients in the national PHRC IC2007-09 study, and this study was authorized by the ethics committees “Comité de Protection des Personnes Sud-Est IV,” references L07–95/L12–171, and “Comité de Protection des Personnes Ile de France,” reference 0811728.
Sample collection and processing
For each case, DNA was extracted from a single tumor sample obtained at diagnosis according to standard procedures. Plasma samples were obtained at diagnosis or during treatment and follow-up (1–6 samples/patient; Supplementary Table S2) by blood sampling directly on standard EDTA tubes and prepared by centrifugation at 2,000 rpm for 10 minutes followed by careful aliquoting and freezing at −80°C within 1 to 24 hours after collection.
cfDNA purification and quantification
cfDNA was extracted from 100 μL to 3.3 mL of plasma using QIAamp Circulating Nucleic Acid Kit (Qiagen) with the Qiavac24s system, according to the manufacturers' recommendation. cfDNA concentration was measured by Qubit fluorometric assay (Invitrogen) with dsDNA HS (High Sensitivity) Assay Kit. The total cfDNA concentration per mL of plasma was calculated. Its quality was defined by bioanalyzer agilent 2100 (Agilent) analysis using the High Sensitivity DNA chip with cfDNA quality expressed as the 200-bp fragment fraction. For 11 plasma samples obtained at diagnosis, cfDNA copy number analysis using Oncoscan has been reported previously (10).
WGS/WES of primary tumors
Primary tumors and paired germline genomic DNA were whole-genome/exome sequenced. For 4 patients, WGS of germline, tumor at diagnosis, and tumor at relapse DNA using an Illumina HiSeq 2500, with 100-bp paired-end reads has been reported previously (11). For the other patients, WES was performed following either an Agilent SureSelect Human All Exon v5 or a Roche Nimblegen SeqCap EZ Exome V3.
Library construction and exome capture of cfDNA
cfDNA libraries were constructed without fragmentation (to account for the mean cfDNA fragment size of 160 bp) using Kapa Library Preparation Kit Illumina platforms (Kapa Biosystems) with Indexed Adapters included in SeqCap EZ Human Exome Kit v3.0 (Nimblegen Roche Sequencing). The manufacturer's protocol was modified with a ligation of 16 hours at 20°C using an adapter:insert molar ratio of 10:1; 10 μL of Pre–LM-PCR oligo 1&2 (5 μmol/L) on 9 cycles of precaptured LM-PCR and 99 μL of SeqCap EZ Purification Beads (1.8×) were used for clean-up of the amplified sample library. Library quantification and quality were determined as above.
For exome capture, SeqCap EZ Exome Enrichment Kit v3.0 (Nimblegen Roche Sequencing) was used according to the manufacturer's protocol. WES was performed on germline, tumor, and plasma samples from 19 patients using Illumina Hi-seq2500 leading to paired-ends 100 × 100 bp. Eight samples were multiplexed for the exome capture and put on 3 lines of high-output flow cell (expected coverage:100×).
For custom deep coverage targeted capture sequencing, a panel was designed encompassing 1,211 SNVs, corresponding to all SNVs observed in WES sequencing. Libraries of cfDNA were constructed using a double capture procedure. Target sequencing was performed on all plasma samples and tumor samples whenever available (PE 150 × 150 bp). Twelve samples were multiplexed for the capture and sequenced with MiSeq V3 reagents (expected coverage: 1,000×).
Bioinformatics pipeline
Following alignment with Bowtie2 (22) allowing up to 4% of mismatches, bam files were cleaned according to the Genome Analysis Toolkit recommendations (23). Targeted bam files were cleaned up without removing duplicates.
WES variant calling was performed using 3 variant callers: GenomeAnalysisTK-3.5 UnifiedGenotyper, HaplotypeCaller, and Samtools-0.1.18 (24). Annovar-v2013-07-29 with cosmic-v64, dbsnp-v137, and RefSeq were used for annotations, and functional prediction was performed using Polyphen2, LRT, MutationTaster, and Sift (25–28). SNVs with a quality <30, a depth of coverage <20 in tumor or plasma samples, or <2 reads supporting the variant were filtered out. Only SNVs within exons of coding genes or splice sites were kept. Then, variants reported in more than 1% of the population in the 1000 genomes (1000gAprl_2012) or Exome Sequencing Project (ESP6500) were discarded to filter out polymorphisms. Finally, synonymous variants were filtered out except those with a COSMIC ID. Variant calling comparison between germline and somatic DNAs allowed us to focus only on tumor-specific SNVs. Only nongermline SNVs supported by >2 reads and a position coverage >20× were taken into account as tumor-specific alterations.
CNAs were analyzed using VarScan.v2.3.5 (29) and DNAcopy-1.42.0 (30) or HMMCopy (31). Estimation of tumor cellularity in analyzed samples (primary tumor or cfDNA) was performed using Sequenza (32).
For the definition of clones, mutated allele fractions (MAF) with similar values were taken into account, while correcting for the copy number status, all similar MAFs (±5%) with values corrected for copy number status defining a clone, as described previously (33). Clonal evolution graphs were generated with Fishplot-0.3 (34). Serial dilutions copy number profiles were generated using CNVkit, and the visualization of log ratio segments was done using matplotlib plotting library (35, 36).
Variants were considered confirmed if observed by two independent sequencing techniques of WES and deep coverage targeted sequencing.
The sequencing data have been deposited at the European Genome-phenome Archive (EGA) under study accession number: EGAS00001002705.
Results
Optimization of library construction for WES cfDNA analysis and bioinformatics pipeline
We have developed an experimental and bioinformatics pipeline to enable WES of cfDNA. Library construction was performed with low cfDNA input (mean, 63 ng; range, 7–100 ng). An additional step of amplification was omitted to reduce false positive rates and further improved quality of sequencing data (37, 38). The step of dual size selection only retained the cfDNA and excluded genomic DNA in case of lymphocyte degradation (Supplementary Fig. S2). These libraries could then be subjected directly to whole-exome or target capture experiments. This optimization enabled a good rate of on-target reads (mean, 82%; range, 69%–88%) and a low percentage of duplicates (mean, 32%; range, 13%–63%) for a WES at 100× coverage (Supplementary Fig. S3), thus enabling a detection of MAF at 5%. The bioinformatics analytic process was established to filter out all germline events based on comparison of primary tumor and cfDNA sequencing results with the germline sequencing results, retaining only somatic alterations for further analysis.
To estimate the sensitivity of cfDNA WES and deep coverage target sequencing, patient cfDNA from one case (Case 14) harboring both significant CNAs and SNVs, and containing a high amount of ctDNA (>90%) was serially diluted with cfDNA from a healthy donor, and dilutions were submitted to both WES and deep coverage capture sequencing. No significant difference of MAFs detected by the two independent techniques for a given SNV was observed, confirming these SNVs and ruling out sequencing errors. Application of the bioinformatics pipeline enabled to establish a threshold for SNV detection of 5% MAFs by WES and 1% by deep coverage target sequencing of cfDNA, below which SNVs were not distinguishable from the background (Supplementary Fig. S4). Furthermore, apart from MYCN amplification, a threshold for detection of CNAs was determined around a proportion of 20% of patient cfDNA in the healthy donor cfDNA (Supplementary Fig. S4).
Comparison between primary tumor and cfDNA at diagnosis
For 19 neuroblastoma patients, WES or WGS data of the primary tumor obtained at diagnosis were available. These data were compared with WES cfDNA analysis at the time of diagnosis.
Analysis based on the Sequenza tool (32) revealed a mean of 73% of tumor cells in the primary tumor sample (range, 15%–98%), and a mean of 60% of ctDNA in the cfDNA fraction (range, 3%–99%; Supplementary Fig. S5). Calling of SNVs could be performed on all diagnostic cfDNA samples, whereas CNAs could not be determined in 1 case (Case 12) due to a ctDNA fraction below the threshold enabling detection of copy number changes.
Following filtering on germline DNA to filter out any constitutional alterations or germline polymorphisms, a total of 861 somatic SNVs were detected by WES analysis in cfDNA, with 353 common to the primary tumor and cfDNA, 97 specific to the primary tumor, and 411 specific to cfDNA (Supplementary Table S3). Overall, a mean of 19 SNVs (range, 9–69) per case common to primary neuroblastoma and cfDNA was observed (Fig. 2A). SNVs observed recurrently at diagnosis in both primary tumor and cfDNA by WES concerned among others the genes ALK (2 cases; Fig. 3D), genes of the MAPK pathway (3 cases), ARHGAP genes (4 cases), NOTCH (1 case), MLL genes (4 cases). SNVs of genes described to play a role in neuroblastoma oncogenesis or other tumors are highlighted in Table 1, and all detected SNVs are listed in Supplementary Table S3. No correlation of MAFs was observed between the solid tumor and cfDNA (R² = 0.017; Supplementary Fig. S6A). A mean of 6 SNVs per case were specific to primary neuroblastoma and 22 specific to the cfDNA, respectively, including among others genes of the MAPK pathway.
Furthermore, of a total of 162 CNAs, 151 CNAs were common to analyses of primary neuroblastoma and cfDNA (Table 1; Figs. 2B and 3A), with only 3 and 8 specific to either. Concerning the known recurrent CNA of neuroblastoma, MYCN amplification, 1p deletion, 11q deletion, and 17q gain were observed in both primary neuroblastoma and cfDNA in 10, 12, 6, and 12 of 19 cases, respectively, with concordant breakpoints. Alterations seen only in cfDNA but not in the primary were 1p deletion and 17q gain in 1 case each. No discrepancies were observed between CNA determined in cfDNA by WES as compared with an analysis using OncoScan for 11 cases reported previously (10).
Amplification of MDM2 and PTPRB are observed for two cases (Case 1 and 17) in primary neuroblastoma and cfDNA. One case (Case 1) demonstrated a focal gain of chromosome 15q harboring the IGF1R gene detected only in cfDNA (Fig. 3B), as confirmed by two independent technics (OncoScan and qPCR; ref. 10).
Analysis of a second time point cfDNA
cfDNA WES was then performed at a second time point (Fig. 1; Supplementary Table S2). For 8 patients, this corresponded to complete or partial remission (CR/PR), whereas for 9 other patients, this corresponded to disease progression or relapse (PD) as documented by clinical and/or radiological evaluations. At the second time point, a significantly higher fraction of ctDNA in the overall cfDNA was observed in patients with PD: for CR/PR patients (8 cases), the mean fraction of ctDNA in the overall cfDNA fraction was 33% (range, 0%–40%), versus a mean of 67% (range, 43%–89%) for PD patients (Supplementary Fig. S5; P = 0.002898). Calling of SNVs could be performed on all second time point cfDNA samples, whereas CNAs could only be determined in 5 cases.
For patients in CR/PR (8 cases), a total of 64 SNVs were detected at the second time point. Among these, a total of 45 SNVs, with a mean of 6 SNVs (range, 2–12) per case, were common to primary neuroblastoma and the diagnostic ctDNA time point (Fig. 2A). For patients in CR/PR at the second time point, the majority of SNVs disappeared, including an ALK mutation (Case 4; Fig. 3D): Of a total of 323 SNVs seen in cfDNA at diagnosis, only 61 SNVs were seen at the second time point. Only in 3 cases was one new SNV observed at this second ctDNA time point (ZNF814, ILDR2, RREB1).
For patients in PD (9 cases), a total of 392 SNVs were detected at the second time point. Among these, 182 were common to SNVs seen in cfDNA at diagnosis, including a ALK mutation (Case 11). A total of 90 SNVs common to both the primary neuroblastoma and the diagnostic ctDNA time point were detected, with a mean of 10 SNVs per case (range, 3–21; Fig. 2A). A total of 198 new SNVs, with a mean of 22 new SNVs per case (range, 0–55) were observed at this second cfDNA time point. Among the new, seemingly relapse/progression–specific SNVs, several SNVs shown previously to play a role in cancer were detected targeting genes such as KRAS (Case 11 relapse; c.G38A; p.G13D; Fig. 3E), ADCY7 (Case 13; c.T2834C;p.V945A), or NCOR1 (Case 15; c.C4367T;p.A1456V). Pathways targeted by these new seemingly relapse-specific SNVs concerned genes of the MAPK pathways and targeted protein kinase A signaling [Ingenuity Pathways Analysis (IPA)]. In addition, in 5 interpretable cases, a total of 16 new CNAs were observed (Fig. 2B), including CDK4 amplification (Case 15, Fig. 3C) or CDK6 amplification associated with BCL11A deletion (Case 16).
Deep coverage capture sequencing of cfDNA at intermediate time points
A panel was then designed encompassing all SNVs observed in any of the primary tumors or cfDNA samples studied by WES, and also taking into account SNVs observed in a relapse tumor for the 4 cases where sequencing data from the relapse tumor were available (n = 1,211; ref. 11), to enable deep coverage capture sequencing of all intermediate time points. In cases with confirmed SNVs, for a given SNV, MAFs observed either by cfDNA WES or by deep coverage capture sequencing did not show a statistically significant difference (P = 0.897, t test) and good positive correlation (R2 = 0.899; Supplementary Fig. S6B), indicating the concordance of the two technics.
Capture sequencing enabled to analyze cfDNA samples obtained at diagnosis and at the second time point for all patients, as well as intermediate time points in 9 patients (Fig. 1) aiming at a high coverage to search for subclonal events. Among the SNVs seen at relapse but not at diagnosis, targeted high coverage sequencing was employed to determine whether they were present in a subclone at diagnosis or not. Interestingly for 17 seemingly relapse-specific SNVs detected by cfDNA WES at relapse but not by tumor or cfDNA WES at diagnosis, higher coverage target analysis of primary tumor or cfDNA from diagnosis showed the presence of these alterations in a minor subclone (MAF mean, 0.84%; range, 0.33%–2.91%), indicating that these minor subclonal populations at diagnosis emerged to major clones during tumor progression and suggesting clonal evolution in 5 of 9 cases with PD (Supplementary Fig. S7). Genes targeted by the alterations in the emerging clones included 10 genes involved in neuritogenesis and 5 in cell cycle and connective tissue development (IPA network analysis; Supplementary Table S1).
Modelization of clonal evolution based on MAFs
Deep coverage capture sequencing also enabled to analyze intermediate time points in 9 patients (Cases 1–9; Supplementary Table S2), aiming at a high coverage to search for subclonal events during treatment or follow-up (Fig. 1). For these 9 cases, the MAFs of SNVs detected in cfDNA were followed in a series of 2 to 6 intermediate plasma samples. For each case, the MAFs corresponded to the clinical disease status. For one case, a decrease of MAFs concurred with disease remission (Case 9; Fig. 4A). For another case (Case 6; Fig. 4B), a decrease of the MAFs was observed during first-line chemotherapy, but an increase of the MAFs was observed corresponding to disease progression (Supplementary Fig. S8A).
Among these 9 cases, the MAFs of all SNVs observed at all time points were used to develop more detailed models of clonal evolution in neuroblastoma in patients with PD (2 patients). Clones were defined on the basis of MAFs with values corrected for copy number status, and the evolution of clones was inferred from MAFs evolving in a similar fashion. This enabled to emit hypotheses with regards to the clonal composition of the tumor at different time points.
In one case with intermediate metastatic progression, then stable disease, a persisting, resistant clone with concomitant disappearance of other clones, was identified by a mutation in HERC2, a gene coding for a ubiquitin protein ligase (Case 1; Fig. 5A). In another case, two minor clones previously detected at the time of diagnosis emerged with an increase of MAFs from 10% to 20% and 1% to 5%, involving the genes ADRM1, BMPR2, CELSR1 (Case 6, Figs. 4B and 5B; Supplementary Fig. S8B).
In conclusion, this novel technology enabled to follow MAFs in sequential cfDNA samples from neuroblastoma patients and indicated different mutational evolutionary patterns. In patients with response to treatment, as expected, evidence of tumor-specific genetic alterations disappeared in cfDNA. On the other hand, in case of persistent disease or relapse, tumor cell–specific alterations could be detected and intermediate analyses revealed either resistance of preexisting or emergence of new clones.
Discussion
Recently, cfDNA has emerged for the analysis of ctDNA, as this surrogate marker can detect tumor cell–specific alterations in a noninvasive setting. The majority of ctDNA studies are based on tools such as digital droplet PCR or targeted sequencing, requiring prior characterization of the biomarkers to be studied (2, 3, 15). However, clonal evolution has been shown to play an important role in many solid cancers, and new clones will not be detected when basing analysis on a comparison with the primary tumors (12, 39, 40), suggesting that WES will be more adapted in particular for cancers with few recurrent biomarkers such as neuroblastoma, despite its lower sensitivity.
Neuroblastoma has been shown to shed important amounts of ctDNA in particular in metastatic cases and high-risk disease and thus represents a relevant model for further ctDNA analysis. In this study, the ctDNA fraction of the overall cfDNA could be calculated by Sequenza (32), confirming a high proportion of ctDNA in the cfDNA fraction (mean, 60%; range, 3%–99%; Supplementary Table S1), whereas for healthy adults, in general, cfDNA quantities of approximately 10 ng/mL of plasma are observed; in high-risk neuroblastoma patients at diagnosis, high amounts of cfDNA (mean, 1,034 ng/mL of plasma; range, 26–13,533) were documented, further confirming the observation of an important shredding of ctDNA into the bloodstream by high-risk neuroblastoma at diagnosis (10, 16). Further studies on larger series will enable to stablish whether the amount of ctDNA is linked to the overall disease burden, the degree of necrosis, or both. Importantly, an increase of the ctDNA fraction in the overall cfDNA was observed at the time of disease progression, versus a decrease with tumor response, suggesting that the ctDNA fraction of cfDNA in general is predictive of the underlying disease status.
WES analysis could be applied successfully to cfDNA samples, enabling calling of SNVs in 36 of 36 and calling of CNAs in 23 of 36 samples. For a given SNV, there was no correlation between the observed MAF in the primary neuroblastoma versus the plasma obtained at the same time point. This might be explained by differences of tumor content between the tumor sample and the ctDNA fraction in the total cfDNA or by spatial heterogeneity with specific SNVs detected in each sample on the other hand. A wide range of MAFs was observed in a given sample, which can be explained by subclonal events on the one hand, or the possibility of differential shredding of ctDNA from different cells.
The WES analysis of primary tumor and cfDNA identified an overlap of only 41% of SNVs at diagnosis, whereas an overlap of 93% was observed for CNAs. This suggests a greater intratumor variability for SNVs and highlights the potential of cfDNA analysis when exploration of the whole spectrum of tumor-specific alterations is of importance, for instance, in the setting of precision medicine. Multiple samples from distinct tumor sites, or tumor and metastatic sites, might give further insight into such heterogeneity. In case of primary tumor-specific alterations, these alterations might correspond to clones that might release less ctDNA, which according to one hypothesis could correspond to less aggressive cells. On the other hand, ctDNA-specific alterations might correspond to more aggressive clones or to those originating from metastatic sites (20, 41, 42). The overall ctDNA MAFs corresponded to the clinical disease status, indicating that ctDNA SNVs might serve as a surrogate marker for disease burden and/or minimal residual disease (43, 44). Furthermore, the emergence of new SNVs coincided with disease progression. As plasma samples in this study were not obtained within a specific systematic surveillance program, we cannot conclude whether cfDNA WES might have served as an early marker for progression. This important question will be addressed in prospective studies.
The study of SNVs and comparison of SNVs between diagnosis and relapse allows the identification of pathways involved in neuroblastoma progression, including MAPK genes, as described previously, and with one pathway specific to relapse involving protein kinase A signaling (45). Furthermore, SNVs that emerged from a subclonal level at diagnosis to a clonal level at relapse preferentially affected pathways of neuritogenesis (Supplementary Table S1; ref. 46).
To study whether seemingly relapse-specific alterations were already present at the time of diagnosis, but in a minor subclone not detected by standard coverage WES, the depth of coverage was increased. This enabled the detection of 17 SNVs, which were initially deemed relapse-specific but which are already present in a minor subclone at diagnosis, underlining the importance of high depth of coverage when searching for subclones. Aiming for a depth of coverage of 1,000×, MAFs >1% were retained as significant. However, increasing the depth of coverage in the current experimental procedure was also associated with PCR duplicates. These duplicates were filtered out during the bioinformatics analysis pipeline, but it cannot be excluded that some identical reads not corresponding to PCR duplicates are also falsely removed. Thus, this protocol could be further improved using a UMI (unique molecular identifiers) prior to library construction to separate real PCR duplicates from alignment duplicates (47).
Altogether, the study of SNVs and analysis of MAFs of a given SNV over time enables the study of clonal evolution and to develop different models of tumor evolution (Figs. 4 and 5; ref. 33). In one model, among different clones derived from the primary clone, one subclone at diagnosis resists to chemotherapy and takes precedence, becoming a major clone at relapse. Indeed, a clone resistant to treatment was identified by the presence of a mutation in a ubiquitin ligase gene. In another model, subclones derived from the primary clone, which remains major, become resistant and collaborate with the primary clone for a relapse. Indeed, it can be hypothesized that these emerging clones might contribute to tumor progression by collaboration with other nonemerging clones (48).
In conclusion, the study of clonal evolution based on cfDNA WES highlights mechanisms of treatment resistance and of escape. These techniques are of high clinical potential in the context of precision medicine and sequential analysis of molecular targets, enabling adaptation of treatment strategies depending on clonal composition. Indeed, biomarker-based clinical trials take into account a molecular profile obtained at a given time point. The evolution of predictive biomarkers over time, which might occur under selective pressure, such as targeted therapies, might now be further evaluated on the basis of noninvasive ctDNA analysis and thus be useful for identification of drug-resistant clones and will in the future lead to further sequential adaptation of biomarker-based treatment approaches.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Authors' Contributions
Conception and design: M. Chicard, O. Delattre, G. Schleiermacher
Development of methodology: M. Chicard, G. Schleiermacher
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): N. Clement, S. Baulande, G. Pierron, E. Lapouble, M. Peuchmaur, N. Corradini, J. Michon, V. Combaret, O. Delattre, G. Schleiermacher
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M. Chicard, N. Clement, S. Baulande, P. Deveau, G. Pierron, E. Lapouble, M. Peuchmaur, N. Corradini, J. Michon, V. Combaret, O. Delattre, G. Schleiermacher
Writing, review, and/or revision of the manuscript: M. Chicard, L. Colmet-Daage, A. Danzon, V. Bernard, S. Baulande, G. Pierron, I. Janoueix-Lerosey, M. Peuchmaur, N. Corradini, A.S. Defachelles, D. Valteau-Couanet, V. Combaret, O. Delattre, G. Schleiermacher
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M. Chicard, N. Clement, M. Bohec, A. Bellini, G. Pierron, I. Janoueix-Lerosey, M. Peuchmaur, A.S. Defachelles, G. Schleiermacher
Study supervision: M. Chicard, G. Schleiermacher
Acknowledgments
This study was supported by the Annenberg Foundation, the Nelia and Amadeo Barletta Foundation (FNAB), and the Association Hubert Gouin Enfance et Cancer. This study was also funded by the Associations Enfants et Santé, Les Bagouz à Manon, Les amis de Claire. Funding was also obtained from SiRIC/INCa (grant INCa-DGOS-4654), from the CEST of Institute Curie, and PHRC IC2007-09 grant. High-throughput sequencing has been performed by the ICGex NGS platform of the Institut Curie supported by the grants ANR-10-EQPX-03 (Equipex) and ANR-10-INBS-09-08 (France Génomique Consortium) from the Agence Nationale de la Recherche ("Investissements d'Avenir" program), by the Canceropole Ile-de-France and by the SiRIC-Curie program - SiRIC Grant “INCa-DGOS- 4654”. We would like to thank Roche Diagnostic France and Magdalena Nawara for their help in protocol modification and the Circulating Biomarkers Team of Institut Curie (Aurore Rampanou, Caroline Garnier-Hego, Jordan Madic, Charles Decraene, and Charlotte Proudhon) for use of their laboratory and their advice.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.