Twenty FL cases diagnosed and clinically followed up in Dokuz Eylül University hospital were included in this study. The patient samples used in this study were approved by the institutional review board of the Medical School at Dokuz Eylül University (protocol no: 2233-GOA). All patients consent to participate in the study. All FL cases were evaluated according to the WHO 2016 classification, using morphological criteria and an appropriate immunohistochemical panel. All of the samples had a follicular pattern of neoplastic lymphoid infiltrate and the follicular dendritic cell meshworks were present in follicles. Based on these findings transformation to DLBCL was excluded. All available clinical and pathological data of the FL patients were obtained from the online clinical database of DEU hospital. Table 1 shows the clinicopathological characteristics (e.g. type of therapy, cancer stage, the percentage of Ki-67⁺ cells) of the FL cases. In FDG PET/CT imaging, malignant tissues can be distinguished with their higher uptake of FDG than normal. Briefly, FDG PET/CT was performed as follows: Six hours after starvation, F-18 FDG was applied intravenously to FL patients. With the purpose of attenuation correction, low dose CT and then PET images were obtained 60 mins post-injection. The number of anatomical regions potentially having tumor involvement was identified from the FDG PET/CT reports based on their SUVmax measurements.

Peripheral blood was collected from newly diagnosed 13 symptomatic and 7 asymptomatic FL cases before treatment. From 12 of 13 symptomatic FL cases, peripheral blood was collected second time after treatment as one FL patient passed away. The treatment and timeframe for peripheral blood samplings are shown in Supplementary Table S1 . For each sample, 12 ml of peripheral blood was withdrawn into K2 EDTA tubes (BD Vacutainer, USA) in the department of hematology at DEU, and placed immediately on ice. Blood samples were transferred into 15 ml tubes, and centrifuged at 1600g for 10 min at 4°C. Supernatants were transferred to 1.5 ml tubes with each containing 1 ml plasma. The tubes were spun down to remove cellular debris at 16.000g for 5 min at 4°C. The supernatants were then transferred into new tubes, and stored in the -80°C freezer until cfDNA isolations. Circulating cell-free DNA (cfDNA) isolation was performed using QIAamp Circulating Nucleic Acid Kit (Qiagen Inc., Germany) according to the manufacturer’s recommendations. The size distribution of plasma cfDNA samples were evaluated with Agilent Bioanalyzer 2100 system to ensure fragment size distribution at around 150 bp ( Supplementary Figure S1 ), which is the typical size for circulating cell-free DNA (22). Granulocytes was obtained from the pellet part of the blood sample by using Ficoll^® Paque Plus (MilliporeSigma, Burlington, MA) per the manufacturer’s instructions. After that, genomic DNA was isolated from granulocytes using PAXgene Blood DNA Kit (Qiagen Inc., Germany).

Concentrations of cfDNA isolated with QIAamp Circulating Nucleic Acid Kit were determined by using Qubit^® DNA Assay Kit and Qubit^® 2.0 Fluorometer. Plasma cfDNA concentrations per 1 ml of plasma were calculated using the following formula: (Qubit concentration measurement as ng/µl * elution volume as µl)/(plasma volume as ml used during cfDNA isolations).

Five µm-thick sections were prepared from tumor tissue blocks that were fixed with formalin and embedded in paraffin for each FL patient. The sections were then placed into 1.5 ml tubes, and stored at 4°C until isolations to prevent nucleic acid degradation (23). Tumor tissue genomic DNA samples were isolated using QIAamp DNA FFPE Tissue Kit (Qiagen Inc., Hilden, Germany) according to the manufacturer’s instructions.

The custom-designed NGS panel including FL and transformed FL-associated 110 genes with somatic mutations was constructed based on previously published studies (24–31). In this targeted NGS panel 575.481 bp genomic region, including the coding sequences and 50 bp flanking intronic sequences of 110 FL-related genes was targeted. FL-associated genes included in the targeted NGS panel are shown in Table 2 .

The targeted sequencing experiments on plasma cfDNA, tumor tissue DNA, and granulocyte gDNA samples from 20 diagnostic FL cases, and 12 plasma cfDNA samples of immunochemotherapy treated FL cases were performed in Novogene Genome Sequencing Company (Hong Kong, China). These target regions were captured with 6385 probes; and Agilent SureselectXT Custom Kit was used for library preparations. Captured regions were paired-end sequenced using a HiSeq system, and 150 bp NGS reads were generated. Targeted genomic regions were covered with > 1500X in-depth for the identification of somatic variants with low variant allele fractions (VAF).

Somatic variants were identified by using a computational bioinformatics pipeline that includes the MuTect2 variant caller ( Supplementary Figure S2 ). The workflow of this pipeline is briefly as follows: First, the quality of raw NGS reads were evaluated with the FASTQC program (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Secondly, the quality of raw NGS reads were improved by applying the AfterQC program (32). Processed NGS reads were then aligned to the human reference genome (hg19) with the BWA-MEM tool (33). Picard tool (https://broadinstitute.github.io/picard/) was used to mark and filter off the PCR duplicates in NGS reads processed with SAMtools (34). MuTect2 variant caller, which can identify somatic variants with a low percentage of VAF, was used to identify somatic variants in plasma cfDNA and tumor tissue DNA samples (35). Patient-matched granulocyte gDNA samples were used to filter off the germline SNPs. Identified variants were annotated to genomic regions with Annovar (36). In the final stage, the aligned NGS reads corresponding to the genomic locations of the identified variants were visually investigated with Integrative Genomics Viewer (IGV) to verify the presence and somatic nature of the mutations (37). During IGV visualizations, the following inclusion criteria were applied for the identification of a confident list of somatic variants: 1) Variants present in the COSMIC database (release v89), and associated with hematopoietic and lymphoid tissues; 2) Missense, nonsense, indel, or splice site mutations; 3) Variants with a VAF > 0.5%; 4) Variants in cfDNA or tumor tissue DNA with a VAF > 4 fold compared to that in the granulocyte gDNA; and 5) Variants with >5 forward and >5 reverse NGS reads.

PCR-Sanger was applied on FL cfDNA and/or tumor tissue DNA to cross-validate selected somatic variants with a percentage of VAF greater than 20%. The DNA nucleotide sequences nearby the selected variants were obtained using the “Get DNA” bioinformatics tool available in UCSC Genome Browser (http://genome.ucsc.edu/). The primer pairs to be used for amplification of the genomic regions including the evaluated variants were designed with the PrimerQuest software (IDT DNA technologies, Coralville, IA). Designed primers were produced by Sentebiolab Biotech (Ankara, Turkey). Optimum annealing temperatures of the primer pairs were determined by gradient PCR. FL cfDNA and tumor tissue DNA samples were amplified with PCR and then send to Eurofins Genomics (Ebersberg, Germany), where PCR purifications and Sanger sequencing reactions were performed with forward and reverse primers. Vector NTI 10.3.0 (Invitrogen, Carlsbad, CA) was used to detect mutations in sequenced FL DNA samples by comparing to the human reference gene sequences. The primers used in PCR-Sanger experiments are shown in Supplementary Table S2 .

RStudio software (https://www.rstudio.com/) was used for linear regression, survival, and signaling pathway analyses. Survival (38) and Survminer (39) R packages were used during survival analyses and generation of Kaplan-Meier graphics. The statistical significance of the differences in survival of FL patients were evaluated with the log-rank test. Cutoff level for high and low values for each evaluated quantitative variable (i.e. cfDNA concentrations, serum LDH levels, Ki-67 index, and FDG PET/CT positive regions) were determined with maximally selected rank statistics using the “surv_cutpoint” function of the Survminer R package. The cutoff levels classifying low and high cfDNA groups were 17 ng/µl for overall survival and 16 ng/µl for progression-free survival analyses. These cutoff levels were 35% for Ki-67 index, 262.0 units per liter for LDH level, and 10 regions for FDG PET/CT positive regions. The R-value of linear regression analyses for the percentage of VAFs was calculated on the Spearman correlation model; and the p-value was generated with t-test. Hypergeometric distribution test was used to determine the statistical significance of enriched pathways by using the ReactomePA package (40). p<0.05 was considered statistically significant for the analyses performed with the RStudio software.

The overall plan of the study involves the following steps: 1) Construction of a custom-designed targeted NGS panel based on previously mutated genes in FL and tFL; 2) Collection of peripheral blood and tumor tissue samples from treatment-naive and immunochemotherapy-treated FL patients; 3) DNA isolations from peripheral blood plasma, granulocyte, and FFPE tumor tissue samples; 4) Identification of cancer-associated mutations through computational bioinformatics analyses; 5) Investigation of the mutual presence of somatic variants in plasma cfDNA and tumor tissue DNA; 6) Survival analyses with identified variants; 7) PCR-Sanger cross-validation of selected variants. The workflow of this study is shown in Supplementary Figure S3 .

Full set of plasma cfDNA, tumor tissue gDNA, granulocyte gDNA from 20 treatment-naive FL cases (n=60) as well as 12 post-treatment plasma cfDNA samples passed the quality parameters. The average basic NGS parameters per FL DNA samples are as follows: Total raw reads: 31.7*10^{^6} ± 12.3*10^{^6}, raw data: 9.5 ± 3.7 gigabyte, the percentage of clean reads to raw reads: 93.9% ± 7.8, the average error rate in all bases: 0.011% ± 0.003, Q20: 96.8% ± 0.6, Q30: 92.5% ± 1.3, and GC content: 48.5% ± 2.2. Basic ultra-deep targeted sequencing statistics are shown in Supplementary Table S3 .

To address whether plasma cfDNA concentrations have any relationship with FL patient prognosis or not, we compared survival of treatment-naive FL patients stratified based on plasma cfDNA concentrations. We observed significantly inferior overall survival ( Figure 1A ) and progression-free survival ( Figure 1B ) of FL patients with high plasma cfDNA concentrations compared to those having low cfDNA concentrations. Importantly, plasma cfDNA concentrations of symptomatic FL patients were significantly higher compared to those of asymptomatic ones ( Figure 1C ). Next, by performing linear regression analyses we evaluated whether plasma cfDNA concentrations correlate with demographic, clinical, or pathological parameters (i.e. age, Ki-67 index, pre-treatment LDH levels, FDG PET/CT positive regions, or SUVmax value) or not. We did not observe any correlation between plasma cfDNA levels and age, Ki-67 index, FDG PET/CT positive regions, SUVmax value, or LDH levels ( Figures 1D–H ).

In 20 treatment-naive FL patients, in total 91 hematological cancer-associated somatic variants were detected in cfDNA and/or tumor tissue DNA samples for 31 of 110 FL-related genes targeted. Among these variants, 30 (33%) of them were common in both plasma cfDNA and tumor tissue DNA, whereas 56 (61.5%) and 5 (5.5%) of them were observed only in tumor tissue DNA or plasma cfDNA, respectively ( Figure 2A ). Of significance, we observed a significant positive correlation (R=0.53, p=0.002) between variant allele fractions of variants in tumor tissues and plasma cfDNA samples ( Figure 2B ). Next, we comparatively analyzed the number and percentage of mutations in cfDNA and/or tumor tissue DNA of symptomatic and asymptomatic FL cases. Among 91 cancer-associated somatic variants, 65 of them were identified in symptomatic, and 26 of them were identified in asymptomatic cases. Twenty seven (41.5%) of these 65 mutations were observed in both cfDNA and tumor tissue DNA, 34 (52.3%) of them were observed only in tumor tissue DNA, and four (6.2%) of them were present only in cfDNA (Figure 2C). In asymptomatic FL patients, 3 (11.6%) of 26 somatic mutations were observed in both cfDNA and tumor tissue DNA, 22 (84.6%) of them only in tumor tissue DNA, and just one (3.8%) of them were present only in cfDNA (Figure 2D). The proportion of somatic mutations present in cfDNA samples was markedly higher in symptomatic cases and most of the variants detected in asymptomatic cases were only present in tumor tissue DNA (Figure 2E). Importantly, we observed at least one variant present in both DNA samples in 10 of 13 (77%) of symptomatic FL cases (Supplementary Figure S4) and in 1 of 7 (14.3%) of asymptomatic cases ( Supplementary Figure S5 ). When we compared the number of mutations per FL case in these two FL groups, we observed that the number of mutations per FL patient was marginally higher in symptomatic cases (average: 5.0; median: 4.0) compared to asymptomatic cases (average: 3.7; median: 3.0) (Figure 2F). In 65 mutations observed in symptomatic cases there were 54 missense, 6 nonsense, 2 indel, and 3 splice site mutations. In 26 mutations observed in asymptomatic cases there were 17 missense, 6 nonsense, 2 indel, and only one splice site mutations ( Figure 2G ). The details of 91 cancer-associated somatic mutations identified through ultra-deep targeted NGS are shown in Supplementary Table S4 .

To address whether genes reported to be recurrently mutated in previous studies are also frequently mutated in our cohort of FL patients, we visually observed the frequency of the mutated genes as well as the presence of these genes in different DNA types by generating waterfall graphics ( Figures 3A, B ). The genes mutated in at least 20% of 20 FL cases are as follows: CREBBP (40%), BCL2 (30%), STAT6 (25%), CARD11 (20%), and EZH2 (20% of cases) ( Figure 3A ). Of note, more than one mutation per gene was observed in some FL cases. CREBBP mutations were generally located on the HAT/KAT domain ( Figure 3C ). Mutations of BCL2 were distributed in different parts of the protein whereas STAT6 mutations were enriched in the STAT binding domain, and the EZH2 mutations were mainly on the SET domain ( Figure 3C ). Among five cancer-associated STAT6 missense variants, three of them (i.e. p.D419G, p.D419N, and p.D523V) were present in both cfDNA and tumor tissue DNA in FL patients. Moreover, all patients with STAT6 mutations were also carrying CREBBP mutations consistent with the observations in a previously reported study (41). In total 6 EZH2 variants were identified in four FL patients, of which two were present in both cfDNA and tumor tissue DNA. In addition, we identified the enriched signaling pathways and biological processes associated with the 31 genes identified to have somatic variations in cfDNA and/or tumor tissue DNA samples. This analysis showed that mutated genes are mainly enriched in “Transcriptional regulation by RUNX1”, “IL-4 and IL-13 signaling”, “Chromatin modifying genes”, and “Signaling by the B-cell receptor” ( Supplementary Figure S6 ). However, it should be noted that the number of genes used as input for the gene ontology analysis is low, which may result in overestimation of the signaling pathways or biological processes identified.

We then asked whether any of the somatic variations identified in diagnostic FL cases persisted after immunochemotherapy treatment in symptomatic FL patients. Sixty-five somatic mutations were identified in cfDNA and/or FFPE tumor tissue DNA in 13 symptomatic FL cases. We could not obtain plasma cfDNA from one of the FL cases (i.e. Case-10) as this patient passed away. Based on targeted ultra-deep NGS analyses, 9 of these 12 symptomatic FL cases had at least one somatic mutation in their plasma cfDNA before treatment started. The total number of somatic mutations detectable in pre-treatment plasma cfDNA of these 9 symptomatic FL cases was 30. We specifically evaluated these 30 mutations in cfDNA samples after immunochemotherapy. In plasma cfDNA samples obtained post-therapy, only one indel (i.e. CARD11 p.S622delS, %VAF: 0.8) among 30 missense, nonsense, and indel variant was detectable based on our stringent variant identification criteria. Of note, there was no significant difference in plasma cfDNA concentrations between treatment-naive and post-therapy FL samples ( Supplementary Figure S7 ).

To address whether the identified mutations predict overall and progression-free survival of FL patients, we applied survival analyses for each of 31 mutated genes with or without mutations in cfDNA and/or tumor tissue DNA samples. When we evaluated mutations detected in cfDNA and/or tumor tissue DNA, mutations in the BCL2 gene predicted a poor overall survival (p = 0.0042) ( Figure 4A ). Similarly, when FL patients are dichotomized based on the presence or absence of BCL2 mutations only in plasma cfDNA, FL patients with BCL2 mutations were significantly (p = 0.0007) associated with inferior overall survival ( Figure 4B ). We also observed poor overall survival (p < 0.0001) of a CCND3-mutated FL patient ( Figure 4C ); however, CCND3 was mutated only in the tumor tissue of this FL case. When FL cases were dichotomized based on the presence of BCL2 mutations in cfDNA and/or tumor tissue DNA, BCL2-mutated FL cases were associated with poor progression-free survival albeit not statistically significant (p = 0.05) ( Figure 4D ). Importantly, FL patients with BCL2 mutations were significantly associated with progression-free survival when only cfDNA mutation status was taken into account ( Figure 4E ). Finally, a CCND3-mutated FL case showed significantly poor progression-free survival ( Figure 4F ). When Fisher’s exact test was used to evaluate the statistical significance of these survival differences, presence of BCL2 mutation in cfDNA was still associated with poor FL survival; whereas, no statistical significance was observed in terms of survival for the CCND3-mutated FL case ( Supplementary Table S5 ). No association with FL patient survival was detected for other mutated genes evaluated. Kaplan-Meier curves for the recurrently (15% or more of FL cases) mutated genes with no association to progression-free survival are shown in Supplementary Figure S8 .

Next we evaluated the relationship between previously established clinicopathological variables (i.e. Ki-67 index, serum LDH levels, PET/CT positivity) routinely used in prognostic evaluations during FL diagnosis on patient survival. As expected, high expression of Ki-67 (p = 0.013), high serum LDH levels (p = 0.00076), and high number of FDG PET/CT positive regions (p = 0.00014) were significantly associated with inferior prognosis ( Supplementary Figure S9 ). To address whether BCL2 genotyping can improve current prognostic evaluations in a non-invasive fashion, we evaluated the co-presence of BCL2 mutations in plasma cfDNA and adverse clinical variables on predicting overall and progression-free survival. Of significance, FL patients carrying BCL2 mutation with high Ki-67 index ( Supplementary Figure S10A ), high LDH level ( Supplementary Figure S10B ), or high number of positive FDG PET/CT regions ( Supplementary Figure S10C ) showed the worst overall survival. Similarly, FL patients with BCL2 mutations and high value for one of these parameters had the most inferior progression-free survival ( Supplementary Figures S10D–F ). The relationship between FL patient survival and different variables (BCL2 mutation status, clinical variables, different treatment types) are altogether shown in Table 3 .

It has been recently shown that mutations in certain genes including those of BTG1 and B2M could identify FL cases with a tendency for early progression and histological transformation (42). In our study, one FL case (Case-16) had the BTG1 p.E46D ( Figure 5A ) and B2M p.Q22* ( Figure 5B ) mutations specifically in plasma cfDNA but not in the tumor tissue DNA of the same FL patient. As expected, patient-matched granulocyte gDNA did not have variant reads, which shows the somatic origin of these mutations. PET/CT imaging showed complete remission 6 months after R-CHOP immunochemotherapy; however, Case-16 showed progressive disease 12 months post-therapy ( Figure 5C ). The observations with this FL case suggest the possibility of identification of high-risk FL patients non-invasively during diagnosis through cfDNA genotyping.

Next, we performed cross-validation of selected variants detected with ultra-deep targeted sequencing using PCR-Sanger. As the threshold for the percentage of VAF for Sanger sequencing is 20%, we chose the variants with a VAF percentage higher than 20%. Cross-validation of selected somatic variants using PCR-Sanger showed 100% consistency with NGS results ( Table 4 ). PCR-Sanger cross-validation of ARID1A p.R1722* variant identified in plasma cfDNA and tumor tissue DNA is shown in Supplementary Figure S11 as a representative.