Longitudinal Circulating Tumor DNA Profiling in Metastatic Colorectal Cancer During Anti-EGFR Therapy

Background Metastatic colorectal cancer (mCRC) is a heterogenous disease with limited precision medicine and targeted therapy options. Monoclonal antibodies against epidermal growth factor receptor (EGFR) have been a crucial treatment option for mCRC. However, proper biomarkers for predicting therapeutic response remain unknown. As a non-invasive test, circulating tumor DNA (ctDNA) is appropriately positioned to reveal tumor heterogeneity and evolution, as it can be used in real-time genomic profiling. To evaluate the significance of ctDNA in monitoring the dynamic therapeutic response and prognosis of mCRC, we detected the baseline and dynamic changes of ctDNA in mCRC patients receiving anti-EGFR therapies. Methods A single-center study was conducted retrospectively. Plasma samples from mCRC patients who received anti-EGFR therapies were collected at baseline and continuous treatment points. The ctDNA was extracted and sequenced with a target panel of tumor-related genes via next-generation sequencing (NGS). Clinical information was also collected and analyzed. Results We conducted dynamic sampling of 22 mCRC patients, analyzed 130 plasma samples, obtained a baseline genomic mutation profile of the patients. In total, 54 variations were detected in 22 plasma samples, with a positive rate of 77.3% (17/22). TP53 was the most mutated gene (59.1%, 13/22), followed by APC (18.2%, 4/22). There was a high concordance rate of genomic characteristics between the tumor tissue test by polymerase chain reaction and ctDNA test by NGS. The mutation discrepancy increased with an extended course of treatment. During remission TP53 and APC were the most frequently decreased clonal mutations and KRAS, NRAS, ERBB2 and PIK3CA were the most decreased subclonal mutations. Both mutation types were increased during progression. The ctDNA decreased earlier than did the responses of computed tomography and traditional tumor markers (carbohydrate antigen 19-9 and carcinoembryonic antigen [CEA]). Lactate dehydrogenase level (P = 0.041), CEA level (P = 0.038), and primary lesion site (P = 0.038) were independent risk factors that influenced overall survival. Moreover, patients with RAS mutations tended to have a worse prognosis (P = 0.072). Conclusions This study demonstrates that ctDNA is a promising biomarker for monitoring the dynamic response to treatment and determining the prognosis of mCRC.

Background: Metastatic colorectal cancer (mCRC) is a heterogenous disease with limited precision medicine and targeted therapy options. Monoclonal antibodies against epidermal growth factor receptor (EGFR) have been a crucial treatment option for mCRC. However, proper biomarkers for predicting therapeutic response remain unknown. As a non-invasive test, circulating tumor DNA (ctDNA) is appropriately positioned to reveal tumor heterogeneity and evolution, as it can be used in real-time genomic profiling. To evaluate the significance of ctDNA in monitoring the dynamic therapeutic response and prognosis of mCRC, we detected the baseline and dynamic changes of ctDNA in mCRC patients receiving anti-EGFR therapies.
Methods: A single-center study was conducted retrospectively. Plasma samples from mCRC patients who received anti-EGFR therapies were collected at baseline and continuous treatment points. The ctDNA was extracted and sequenced with a target panel of tumor-related genes via next-generation sequencing (NGS). Clinical information was also collected and analyzed.
Results: We conducted dynamic sampling of 22 mCRC patients, analyzed 130 plasma samples, obtained a baseline genomic mutation profile of the patients. In total, 54 variations were detected in 22 plasma samples, with a positive rate of 77.3% (17/22). TP53 was the most mutated gene (59.1%, 13/22), followed by APC (18.2%, 4/22). There was a high concordance rate of genomic characteristics between the tumor tissue test by polymerase chain reaction and ctDNA test by NGS. The mutation discrepancy increased with an extended course of treatment. During remission TP53 and APC were the most frequently decreased clonal mutations and KRAS, NRAS, ERBB2 and PIK3CA were the most decreased subclonal mutations. Both mutation types were increased during progression. The ctDNA decreased earlier than did the responses of computed tomography and traditional tumor markers (carbohydrate antigen 19-9 and carcinoembryonic antigen [CEA]). Lactate dehydrogenase level (P = 0.041), CEA level (P = 0.038), and primary lesion site (P = 0.038) were independent risk factors that

INTRODUCTION
Colorectal cancer (CRC) is the third most common cancer and second most frequent cause of cancer-related death worldwide (1). Most CRC patients are diagnosed in an advanced stage; even patients diagnosed in an early stage will develop advanced disease. Palliative chemotherapy has been the mainstay treatment for metastatic CRC (mCRC), and the overall survival (OS) rate of mCRC patients is less than 3 years (2). The emergence and application of targeted therapies have greatly improved OS (3). Cetuximab is a chimeric human/mouse immunoglobulin G1 monoclonal antibody that targets the human epidermal growth factor receptor (EGFR) protein.
Several international multicenter clinical studies have shown that cetuximab extends median survival to approximately 30 months in RAS and BRAF wild-type mCRC cases (4)(5)(6).
It is important to note that 40-60% of mCRC patients with initial wild-type RAS and BRAF genotypes develop drug resistance after prolonged exposure to cetuximab (7)(8)(9). There is evidence that mitogen-activated protein kinase (MAPK) signal transduction pathway-related genes, c-Met gene amplification, and secondary changes of other genes in the human EGF family may be important mechanisms underlying resistance to anti-EGFR monoclonal antibodies (10). Misale et al. (11) showed at both the cellular level and in the clinic that secondary KRAS mutations may be the mechanism responsible for drug resistance following EGFR blockade. Genomic analyses of biopsied tissue after the development of anti-EGFR therapy resistance have shown multiple mutations in KRAS, NRAS, BRAF, and phosphoinositide 3-kinase catalytic subunit alpha (PIK3CA) genes (12,13). In addition, abnormal changes in genes in the HER family contribute to anti-EGFR resistance. Such changes include mutations in the EGFR extracellular domain and HER-2 amplification, which is common in breast and gastric cancers. Therefore, assessing genomic alterations will help identify potential drug resistance, allowing physicians to adjust treatment decisions in a timely manner.
Circulating tumor DNA (ctDNA) originates from the apoptotic and necrotic turnover of cancer cells. Its genomic profile corresponds with the tumor DNA from which it was derived. Previous studies have shown that ctDNA is enriched in the plasma of cancer patients and its characteristics are representative of the entire tumor genome (14)(15)(16). By defining the genomic features in patient plasma with next-generation sequencing (NGS), ctDNA-based liquid biopsy is a convenient, minimally invasive test with reproducible results. It also complements the limitations of tissue evaluation and helps to monitor the molecular changes that occur during cancer evolution (17)(18)(19). In addition, the consistency of the RAS (including KRAS and NRAS) gene between tissue samples and liquid biopsy is approximately 93% (20). Therefore, ctDNA detection has been used to make treatment decisions in various types of cancer such as colorectal, lung, and gastroesophageal cancers (21)(22)(23).
Previous studies have revealed the importance of ctDNA in CRC. For patients who undergo surgery, ctDNA levels have been used to detect minimal residual disease and predict prognosis (24,25). Moreover, the use of ctDNA in dynamic monitoring of the treatment response has been actively explored. There is growing evidence supporting the significance of ctDNA in the therapeutic response and drug resistance. Detectable ctDNA levels at baseline and new emerging ctDNA at follow-up treatment are associated with a poor prognosis (26). Moreover, decreased ctDNA levels reflect sensitivity to anti-EGFR therapies (27,28). However, more studies are needed to identify the ctDNA features and dynamic changes in mCRC patients.
In this study, we identified ctDNA profiling at baseline and dynamic changes during anti-EGFR treatments. We also explored the relationships between ctDNA abundance and clinical characteristics, prognosis, and therapeutic evaluation.

Patients
Patients pathologically diagnosed with mCRC at Fudan University Shanghai Cancer Center (Shanghai, China) from October 12, 2016 to March 20, 2020 were included in this study retrospectively. The inclusion criteria were as follows: diagnosis of mCRC; presence of at least one measurable or unmeasurable but evaluable lesion (described according to Response Evaluation Criteria in Solid Tumors [RECIST] 1.1); presence of polymerase chain reaction (PCR)-confirmed wildtype KRAS (exon 2/3/4), NRAS (exon 2/3/4), and BRAF (exon 15) genotypes in tumor tissue before the receipt of anti-EGFR therapy; no history of severe heart or liver disease, psychiatric disorders, hemorrhage, or perforation of the digestive tract; and an Eastern Cooperative Oncology Group performance status of 0/1 at 3 days before treatment. Exclusion criteria were as follows: presence of mCRC combined with other types of cancer. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The study protocol was approved by the Ethics Committee of Fudan University Shanghai Cancer Center (Shanghai, China). All patients provided written informed consent to participate. Tumor burden was measured to evaluate the clinical response by computed tomography (CT)

Targeted Capture Sequencing
Cell-free DNA libraries were prepared using the KAPA Hyper Prep Kit (KAPA Biosystems Inc., Wilmington, MA, USA) in accordance with the manufacturer's protocol. They were individually barcoded with unique molecular identifiers. In brief, 30-60 ng ctDNA were subjected to end-repairing, Atailing, and ligation with indexed adapters. Then, the libraries were PCR-amplified and purified for target enrichment. The concentration and size distribution of each library were determined using a Qubit 3.0 fluorometer (Thermo Fisher Scientific) and a LabChip GX Touch HT Analyzer (PerkinElmer, Waltham, MA, USA), respectively. For targeted capture, indexed libraries were subjected to probe-based hybridization with a customized NGS panel that included 61 cancer-related genes. The probe baits were used to individually synthesize 5′ biotinylated 120 base pair (bp) DNA oligonucleotides (IDT, Coralville, IA, USA). Repetitive elements were filtered out from intronic baits according to annotations by UCSC Genome RepeatMasker (29). The xGen ® Hybridization and Wash Kit (IDT) was employed for hybridization enrichment. Briefly, 500 ng indexed DNA libraries were pooled to obtain 2 mg DNA. Pooled DNA samples were mixed with Human Cot-1 DNA and xGen Universal Blockers-TS Mix and dried in a SpeedVac system. Hybridization Master Mix was added to each sample. The mixtures were incubated in a thermal cycler at 95°C for 10 min, then combined with 4 mL probes and incubated at 65°C overnight. Target regions were captured in accordance with the manufacturer's instructions. The concentration and fragment size distribution of the final library were determined using a Qubit 3.0 fluorometer (Thermo Fisher Scientific) and LabChip GX Touch HT Analyzer (PerkinElmer), respectively. The captured libraries were loaded onto a NovaSeq 6000 platform (Illumina, San Diego, CA, USA) for 100 bp paired-end sequencing with a mean sequencing depth of 36000.
Raw data were mapped to the reference human genome hg19 using the Burrows-Wheeler Aligner. In-house developed software was used to generate duplex consensus sequences based on dual unique molecular identifiers integrated at the ends of the DNA fragments. To improve specificity, particularly for variants with low allele frequency in the ctDNA, an in-house loci-specific variant detection model based on a binomial test was applied. The variants were subsequently filtered according to their supporting count, strand bias status, base quality, and mapping quality. In addition, variant calling was optimized to detect variants in short tandem repeat regions. Single-nucleotide polymorphisms (SNPs) and indels were annotated by ANNOVAR against the following databases: dbSNP (v138), 1000Genome, and ESP6500 (population frequency > 0.015). Only missense, stop-gain, frameshift, and non-frameshift indel mutations were kept. Copy number variations and gene rearrangements were detected as described previously (30). We calculated the sum of the variant allele frequency (VAF) for each sample. In this manner, the sum of VAF in percentages represented most of the ctDNA detected at each time point. Because there are no established cutoffs for clinically significant treatment-induced changes in ctDNA, we predefined molecular progression as an increase in the mean VAF by at least 25% or new emerging variant allele if the VAF was negative at baseline. We predefined molecular remission as a decrease in the mean VAF by at least 50% for patients in whom the VAF was positive at baseline. A mutation was defined as "subclone" if the VAF was less than 25% of the highest in the sample or as "clone" if the VAF was above this threshold, according to the method used in a previous study (15).

Statistical Analyses
Numerical diversity between subgroups was assessed using the Wilcoxon-Mann-Whitney test. Survival results were assessed by Kaplan-Meier survival analysis paired with the log-rank test and Cox proportional hazards modeling. For all tests, P < 0.05 was considered statistically significant. All statistical analyses were performed using SPSS software (v. 21.0; SPSS Inc., Chicago, IL, USA) or GraphPad Prism (v. 8.0; La Jolla, CA, USA).

Mutation Profiles at Baseline
For the 22 enrolled patients, ctDNA was extracted from 130 plasma samples and sequenced by NGS. The genomic features of 61 genes (Tables S1), including copy number variations and mutations, were evaluated as treatment proceeded. The mutation profiles at baseline are shown in Figure 1. Using targeted capture sequencing, 54 variations were detected in 22 plasma samples, with a positive rate of 77.3% (17/22). Fifteen genes with different types of variation were identified, including missense, frameshift, stop-gain, gain, noncoding, and multiple variations ( Figure 1A). The RAS mutation discrepancy was also compared among treatments ( Figure 1B). For patients who received cetuximab as first-line treatment, the RAS mutation discrepancy was 13.3% (2/15). Both of these patients also had NRAS mutations. The mutation sites were NRAS p.Q61K (0.31%), NRAS p.G13R (0.07%), and NRAS p.G12R (0.37%). For patients who received cetuximab as second-line treatment, the RAS mutation discrepancy was 33.3% (2/6). One patient had a KRAS p.G12V mutation (2.17%) and the other patient had both KRAS p.Q61H (0.02%) and NRAS p.G13C (0.03%) mutations. The only patient who received cetuximab as third-line treatment had a KRAS mutation. The mutation sites included KRAS p.Q61Hc.183A>T (0.05%), KRAS p.Q61Hc.183A>C (0.91%), and KRAS p.G12A (0.58%). The clonal and subclonal landscapes were detected at baseline ( Figure 1C). Subclonal mutations were found in 31.8% (7/22) of the patients. The three most common clonal mutation genes were TP53, APC, and BRAF, while the three most common subclonal mutation genes were PIK3CA, KRAS, and NRAS.
To further evaluate the consistency of the gene mutation profile in ctDNA and clinical parameters, the following is a description of a typical case. This patient (P1) had a primary tumor on the left side with synchronous liver metastasis and tumor resection before receiving cetuximab. Throughout the course of treatment, this patient received three different lines of treatment. For the first line, he received cetuximab in combination with FOLFOX for 6 months and cetuximab in combination with leucovorin and 5-fluorouracil for another 10 months, and then got progressive disease and changed to second line therapy. This patient declined intensive chemotherapy including venous 5-fluorouracil at that time, so he received cetuximab in combination with irinotecan as second-line treatment. After progression at February 2019, he refused any treatment. And then he began to receive oral regorafenib as third line treatment at May 2019. He did not received bevacizumab in his therapeutic process for the financial reason. Figure 3A illustrates the overall treatment procedure in P1 and the corresponding lesions on CT. The changes in tumor diameter and tumor antigen biomarkers (carbohydrate antigen 19-9 [CA19-9] and CEA) are shown in Figure 3B. Figure 3C illustrates the serial ctDNA testing in P1, showing the emergence of clonal alterations through the treatment process. As shown in Figure S2, the ctDNA decreased to a very low level in June 2017, which was 3 months earlier than the responses of CT and traditional tumor markers (CA199 and CEA). The first progressive disease was observed in July 2018. In contrast, increased ctDNA was detected in May 2018, which was 2 months earlier than the CT results. These results support the use of ctDNA for evaluating treatment efficacy in advanced CRC cases.

Association Between Genomic Features of ctDNA and Prognosis
To investigate the significance of ctDNA in prognosis, the associations of genomic features of ctDNA with PFS and OS were analyzed. For PFS, univariate survival analysis revealed that the histological type (P = 0.045) and RAS status at baseline (P = 0.002) were poor prognostic indicators ( Table 2). For OS, univariate survival analysis revealed that LDH level (P = 0.047), CEA level (P = 0.008), metastasis sites involved (P = 0.012), histological type (P = 0.003), and RAS status at baseline (P = 0.002) were poor prognostic indicators ( Table 3). The variable factors, including LDH level, CEA level, metastatic sites involved, histological type, differentiation, mean VAF, primary lesion site, absence or presence of primary tumor resection, and RAS status at baseline, were included in the multivariate analyses. Although no factor was associated with PFS, LDH level (P = 0.041), CEA level (P = 0.038), and primary lesion site (P = 0.038) were independent risk factors. In addition, patients with a RAS mutation tended to have a worse prognosis (P = 0.072).

DISCUSSION
Considering the growing evidence that ctDNA sequencing could represent a valuable resource for genomic discovery, we analyzed the ctDNA profiles of 22 CRC patients from a single-center retrospectively. We found a high similarity of genomic alterations in ctDNA and tumor tissue, similar to that described in a previous report (31). We also identified the baseline characteristics and dynamic changes in ctDNA mutations during anti-EGFR treatment. These results suggest that ctDNA is a stable biomarker available for auxiliary clinical diagnosis, as well as for evaluating CRC tumor progression. Currently, CT and MRI scans are recommended as the main diagnostic and surveillance methods for mCRC patients (32). However, only enlarged tumors can be identified in this manner. Such tumors always exhibit drug resistance, which limits treatment effectiveness (33). Thus, predictive and prognostic markers that represent therapeutic resistance at an early stage are urgently needed. CEA and LDH are often used for auxiliary diagnosis and therapeutic evaluation. Yet, they are insufficient for reflecting genomic alterations. Thus, there is an unmet clinical need for a biomarker that more accurately reflects therapeutic efficacy and dynamic changes during therapy. Liquid biopsy, particularly ctDNA from plasma, has high sensitivity and specificity in early cancer detection (34), so it may be useful for assessing dynamic changes in genes. Furthermore, because tumor tissue sequencing often relies on archival tissue obtained prior to the development of metastatic disease, ctDNA profiling may more readily facilitate the analysis of patients with metastatic disease by better capturing the presence of tumor heterogeneity.
Our study focused mainly on therapeutic resistance. Because tissue-based sequencing compendia depend mainly on treatmentnaïve tumors in the early stages of development, the results generally cannot focus on acquired resistance. Conversely, ctDNA sequencing can more easily provide non-invasive access to patients with advanced tumors and offer unique insights into resistance mechanisms that emerge under the selective pressures of different therapies. Our analysis identified clonal and subclonal gene mutations, the frequencies of which were decreased during remission and increased during progression. These results confirmed that the recurrent alterations of these previously identified gene alterations in ctDNA were decreased during remission. In contrast, the ctDNA appeared again or the corresponding number of gene alterations was increased during progression, which is possibly associated with tumor progression.
Notably, the acquired mutations of KRAS or NRAS can be detected in the ctDNA of mCRC patients who initially exhibited wild-type genotypes and thus received anti-EGFR therapies. The mutation frequency fluctuated in a dynamic manner according to treatment efficacy, consistent with our dynamic monitoring results. After changes to the secondary therapeutic approach without anti-EGFR antibodies, the mutation frequency can greatly decrease and may disappear. This provides an opportunity to re-challenge the tumor with anti-EGFR antibodies. Some retrospective studies have shown that patients with wild-type RAS and BRAF genotypes can benefit such re-challenging (35). Furthermore, the CAVE study  demonstrated that cetuximab-based re-challenge therapy in RAS wild-type mCRC, according to the results of ctDNA analysis, could be used for patient selection and may improve OS (36). We are also performing a prospective phase II study to evaluate the significance of ctDNA for treatment decision-making in patients with mCRC after failed first-line cetuximab treatment (NCT04831528). This study has some limitations. First, the small number of patients included may have diluted the importance of ctDNA as a predictive and prognosis marker in mCRC. According to a previous study, ctDNA has an important role in advanced solid tumors and can be used to predict treatment efficacy in perioperative CRC (15,22). It could also be a significant marker in mCRC. This limitation can be overcome by validation in a larger study. Second, genomic alterations in ctDNA were not detected in 15% of cases, which was similar to the rates of ctDNA detection in other CRC series. Some patients may not have had alterations in genes covered by the NGS assay. However, in most cases, the lack of detected genomic alterations in ctDNA was generally caused by other factors, including low tumor burden, absence of ctDNA shedding by some tumors, and timing of blood collection. Some other approaches, such as multiomics-like methylation, exosomes, circulating microRNA, metabonomics, and/or molecular imaging methods (37,38), could possibly be used to detect ctDNA at lower thresholds with greater accuracy and provide more practical value for ctDNA detection.

CONCLUSION
This study demonstrates that ctDNA may be a reliable biomarker to assist in the prognostic evaluation and assessment of treatment efficacy in advanced CRC patients. Assessments of dynamic changes in ctDNA in mCRC patients can identify baseline values for prognostic evaluation and help with clinical decision-making.

DATA AVAILABILITY STATEMENT
The data presented in the study are deposited in the Biological ProjectLibrary repository, accession number PRJCA008093, https:// ngdc.cncb.ac.cn/bioproject/browse/PRJCA008093.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of Fudan University Shanghai Cancer Center. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
WY contributed to the data collection, data analysis and writing the manuscript. JZ assisted in data analysis and editing this manuscript. YL assisted in data analysis and data collection. RL and XYZ assisted in data collection. ZY and SC were mainly responsible for genetic testing of samples. WG, WL, MH and XDZ offered part of cases. ZC contributed to the research design, sample collection and manuscript revisions. All authors contributed to the article and approved the submitted version.